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Abstract 

We provide ingredients and recipes for computing signals of TeV- 
scale Dark Matter annihilations and decays in the Galaxy and 
beyond. For each DM channel, we present the energy spectra of 
e ± , p, J, 7, v e ^ T at production, computed by high-statistics simula- 
tions. We estimate the Monte Carlo uncertainty by comparing the 
results yielded by the Pythia and Herwig event generators. We 
then provide the propagation functions for charged particles in the 
Galaxy, for several DM distribution profiles and sets of propaga- 
tion parameters. Propagation of e ± is performed with an improved 
semi-analytic method that takes into account position-dependent 
energy losses in the Milky Way. Using such propagation functions, 
we compute the energy spectra of e ± ,p and d at the location of the 
Earth. We then present the gamma ray fluxes, both from prompt 
emission and from Inverse Compton scattering in the galactic halo. 
Finally, we provide the spectra of extragalactic gamma rays. All 
results are available in numerical form and ready to be consumed. 
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1 Introduction 



Cosmology and astrophysics provide several convincing evidences of the existence of Dark 
Matter [1, 2, 3]. The observation that some mass is missing to explain the internal dynamics 
of galaxy clusters and the rotations of galaxies dates back respectively to the '30s and the 
'70s [4]. The observations from weak lensing [5], for instance in the spectacular case of the 
so-called 'bullet cluster' [6], provide evidence that there is mass where nothing is optically 
seen. More generally, global fits to a number of cosmological datasets (Cosmic Microwave 
Background, Large Scale Structure and also Type la Supernovae) allow to determine very 
precisely the amount of DM in the global energy-matter content of the Universe at f2rjM^ 2 = 
0.1123 ± 0.0035 [7] 1 . 

All these signals pertain to the gravitational effects of Dark Matter at the cosmological 
and extragalactical scale. Searches for explicit manifestation of the DM particles that are 
supposed to constitute the halo of our own galaxy (and the large scale structures beyond it) 
have instead so far been giving negative results, but this might be on the point of changing. 

Indirect searches for Dark Matter aim at detecting the signatures of the annihilations 
or decays of DM particles in the fluxes of cosmic rays, intended in a broad sense: charged 
particles (electrons and positrons, protons and antiprotons, deuterium and antideuterium), 
photons (gamma rays, X-rays, synchrotron radiation), neutrinos. Pioneering works have 
explored this as a promising avenue of discovery since the late-70's: gamma rays from 
annihilations were first considered in [8, 9, 10] and then revisited in [11], antiprotons in [12, 
13] and then more systematically in [11, 14, 15], positrons in [12, 11, 14, 16], antideuterons 
have been first discussed in [17, 18, 19], radio-waves from synchrotron radiation from DM 
in [20, 21, 22, 23] and later in [24] (which questions the approach in [22, 23]), extragalactic 
gamma rays have been first discussed in [25] . Inverse Compton gamma rays from DM have 
been only relatively recently considered as a possible signal (see e.g. [26, 27, 28]). In general, 
a key point of all these searches is to look for channels and ranges of energy where it is 
possible to beat the background from ordinary astrophysical processes. This is for instance 
the basic reason why searches for charged particles focus on fluxes of antiparticles (positrons, 
antiprotons, antideuterons), much less abundant in the Universe than the corresponding 
particles. 

A well spread theoretical prejudice wants the DM particles to be thermal relics from the 
Early Universe. They were as abundant as photons in the beginning, being freely created 
and destroyed in pairs when the temperature of the hot plasma was larger then their mass. 
Their relative number density started then being suppressed as annihilations proceeded 
but the temperature dropped below their mass, due to the cooling of the Universe. Finally 
the annihilation processes also froze out as the Universe expanded further. The remaining, 
diluted abundance of stable particles constitutes the DM today. As it turns out, particles 
with weak scale mass (~ 100 GeV — 1 TeV) and weak interactions [1, 2] could play the 
above story remarkably well, and their final abundance would automatically (miracolously?) 
be the observed Odm- While this is not of course the only possibility, the mechanism is 
appealing enough that a several- GeV-to-some- TeV scale DM particle with weak interactions 
(WIMP) is often considered as the most likely DM candidate. 

In any case, this mass range (TeV-ish DM) has the best chances of being thoroughly 
explored in the near future by charged particle and photon observatories, also in combina- 

^^Here 57dm = Pum/Pc is defined as usual as the energy density in Dark Matter with respect to the 
critical energy density of the Universe p c — 3Hq/8ttGn, where Hq is the present Hubble parameter, h is 
its reduced value h = i/o/100 kms _1 Mpc _1 . 
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tion with direct DM searches (aiming at detecting the nuclear recoil produced by a passing 
DM particle in ultra-low background underground detectors) and, possibily, production at 
LHC collider. It is therefore the focus of our attention. 

Supposing (and hoping) therefore that anomalous features are detected in the fluxes of 
cosmic rays, it will be crucial to be able to 'reverse engineer' them to determine which Dark 
Matter is at their origin. Moreover, it will be useful to be able to quickly compute which 
other associated signals are implied by a possible positive detection and have to be looked 
for in other channels. Only via a cross-correlation of multi-messenger signals a putative 
detection of DM will be confirmed or disproved. More generally, in order to compute the 
predicted signatures of a given model of Dark Matter, a number of particle physics and 
astrophysics ingredients are needed. These ingredients are what we aim to provide. 

This work does not contain any new theoretical proposal nor any new study of (existing 
or foreseen) data. It contains instead all the phenomenological ingredients that allow to 
perform the analyses sketched above in the most general possible way. 

More precisely, the rest of this compilation is organized as follows. In Section 2 we start 
by recalling the most commonly used DM distribution profiles in the Milky Way, that we 
will adopt for the computation of all signals. In Section 3 we discuss the production of the 
fluxes of Standard Model particles from DM annihilations (and decays): we compare the 
Pythia and Herwig Monte Carlos and quantify the uncertainties. Section 4 deals with 
the propagation in the Galaxy and the resulting fluxes of charged cosmic rays from DM: 
electrons, positrons, antiprotons, antideuterons. Section 5 deals with the basics of prompt 
gamma rays from Dark Matter annihilations (or decays). Section 6 discusses the 'secondary' 
radiation from e produced by DM annihilations or decays: Inverse Compton Scattering 
(ICS) 7-rays and synchrotron radiation. Section 7 presents the results on extragalactic 
gamma rays. One signature that we do not discuss here is that of neutrino fluxes from the 
annihilation/decays of DM accumulated in the center of the Sun: we refer the reader to [99, 
34], where they have been discussed in a spirit very similar to the one of the present work, 
and we intend to upgrade those former results in the light of the current developements in 
upcoming work. 

Several of these parts contain innovations with respect to the existing literature. For 
instance, the comparison among Monte Carlo generators; the propagation halo functions 
for e , which allow to take into account, in a semi- analytic way, point-dependent energy 
losses and therefore to compute much more precisely the fluxes of charged cosmic rays and 
(above all) ICS photons; the introduction of a formalism in terms of (other) halo functions 
to compute the flux of IC 7 rays; the study of the impact of different choices for the model 
of extragalactic background light on the predicted fluxes of extragalactic gamma rays... 

All our numerical results are available at the website referenced in [29]. So finally in 
Section 8 we give a summary of these provided numerical ingredients and we list the entry 
points in the text for the main recipes. 

Of course, many refined numerical tools which allow to (directly or indirectly) com- 
pute Dark Matter indirect detection signatures have been developed in the latest decades. 
Among them, GALPROP [30], DarkSUSY [31], MicrOMEGAs [32], IsaTools [33], Wimp- 
Sim [34]... Rather than focusing on a particular DM model, we try to be model- independent 
and parameterize the observables in terms of the DM mass, of the DM decay or annihilation 
rates and channels, as well as in terms of a few uncertain astrophysical assumptions. We 
prefer, whenever possible, a semi-analytical treatment that allows us to keep track of the 
approximations and choices that we make along the way. Also, we aim at providing the 
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reader with ready-to-use final products, as opposed to the generating code. We make an 
effort to extend our results to large, multi-TeV DM masses (recently of interest because 
of possible multi-TeV charged cosmic ray anomalies) and small, few-GeV DM masses (re- 
cently discussed because of hints from DM direct detection experiments), at the edge of the 
typical WIMP window. Above all, our aim is to provide a self-consistent, independently 
computed, comprehensive set of results for DM indirect detection. Whenever possible, we 
have compared with existing codes, finding good agreement or improvements. 



2 Dark Matter distribution in the Galaxy 

For the galactic distribution of Dark Matter in the Milky Way we consider several possi- 
bilities. The Navarro, Frenk and White (NFW) [35] profile (peaked as r _1 at the Galactic 
Center (GC)) is a traditional benchmark choice motivated by N-body simulations. The 
Einasto [36, 37] profile (not converging to a power law at the GC and somewhat more 
chubby than NFW at kpc scales) is emerging as a better fit to more recent numerical sim- 
ulations; the shape parameter a varies from simulation to simulation, but 0.17 seem to 
emerge as a central, fiducial value, that we adopt. Cored profiles, such as the truncated 
Isothermal profile [38, 39] or the Burkert profile [40], might be instead more motivated by 
the observations of galactic rotation curves, but seem to run into conflict with the results of 
numerical simulations. On the other hand, profiles steeper that NFW had been previously 
found by Moore and collaborators [41]. 

As long as a convergent determination of the actual DM profile is not reached, it is 
useful to have at disposal the whole range of these possible choices when computing Dark 
Matter signals in the Milky Way. The functional forms of these profiles read: 
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Numerical DM simulations that try to include the effects of the existence of baryons have 
consistently found modified profiles that are steeper in the center with respect to the DM- 
only simulations [42]. Most recently, [43] has found such a trend re-simulating the haloes 
of [36, 37]: steeper Einasto profiles (smaller a) are obtained when baryons are added. 
To account for this possibility we include a modified Einasto profile (that we denote as 
EinastoB, EiB in short in the following) with an a parameter of 0.11. All profiles assume 
spherical symmetry 2 and r is the coordinate centered in the Galactic Center. 

Numerical simulations show that in general halos can deviate from this simplest form, and the isodensity 
surfaces are often better approximated as triaxial ellipsoids instead (e.g. [44]). For the case of the Milky 
Way, however, it is fair to say that at the moment we do not have good observational determinations of its 
shape, despite the efforts already made studying the stellar tidal streams, see [45]. Thus the assumption 
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Figure 1: DM profiles and the corresponding parameters to be plugged in the functional forms 
of eq. (1). The dashed lines represent the smoothed functions adopted for some of the computations 
in Sec. 4-1-3. Notice that we here provide 2 (3) decimal significant digits for the value ofr s (p s ): 
this precision is sufficient for most computations, but more would be needed for specific cases, such 
as to precisely reproduce the J factors (discussed in Sec. 5) for small angular regions around the 
Galactic Center. 

Next, we need to determine the parameters r s (a typical scale radius) and p s (a typical 
scale density) that enter in each of these forms. Instead of taking them from the individual 
simulations, we fix them by imposing that the resulting profiles satisfy the findings of 
astrophysical observations of the Milky Way. Namely, we require: 

- The density of Dark Matter at the location of the Sun r = 8.33 kpc (as determined 
in [48]; see also [49] 3 ) to be p = 0.3 GeV/cm 3 . This is the canonical value routinely 
adopted in the literature (see e.g. [1, 2, 51]), with a typical associated error bar of 
±0.1 GeV/cm 3 and a possible spread up to 0.2 — > 0.8 GeV/cm 3 (sometimes refereed 
to as 'a factor of 2'). Recent computations have found a higher central value and 
possibly a smaller associated error, still subject to debate [52, 53, 54, 55]. 

- The total Dark Matter mass contained in 60 kpc (i.e. a bit larger than the distance to 
the Large Magellanic Cloud, 50 kpc) to be M 60 = 4.7x 1O 11 M . This number is based 
on the recent kinematical surveys of stars in SDSS [56]. We adopt the upper edge of 
their 95% C.L. interval to conservatively take into account that previous studies had 
found somewhat larger values (see e.g. [57, 58]). 

The parameters that we adopt and the profiles are thus given explicitly in fig. 1. Notice that 
they do not differ much (at most 20%) from the parameter often conventionally adopted in 
the literature (see e.g. [2]), so that our results presented below can be quite safely adopted 
for those cases. 

of spherical symmetry, in absence of better determinations, seems to be still well justified. Moreover, it is 
the current standard assumption in the literature and we therefore prefer to stick to it in order to allow 
comparisons. In the future, the proper motion measurements of a huge number of galactic stars by the 
planned GAIA space mission will most probably change the situation and give good constraints on the 
shape of our Galaxy's DM halo, e.g. [46], making it worth to reconsider the assumption. For what concerns 
the impact of non-spherical halos on DM signals, charged particles signals are not expected to be affected, 
as they are sensistive to the local galactic environment. For an early analysis of DM gamma rays al large 
latitudes see [47]. 

3 The commonly adopted value used to be 8.5 kpc on the basis of [50]. 
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As well known, the profiles differ most in the inner part of the galactic halo, close to the 
galactic center, while they are quite self-similar above a few kiloparsecs, and in particular 
around the location of the Earth. As a consequence, DM signals from the inner Galaxy 
(e.g. gamma ray fluxes from regions a few degrees around the GC) will be more sensitive 
to the choice of profile than DM signals that probe the local environment (e.g. the fluxes 
of high energy positrons, produced at most a few kpc away from the Earth) or that probe 
regions distant from the GC (e.g. gamma rays from high latitudes). 

Notice that we will not consider throughout the paper (with the exception of Sec. 7) 
the potential contribution from galactic DM substructures. It is well known that the Cold 
DM paradigm predicts DM halos to contain in a hierarchical fashion a copious number of 
subhalos, something which is clearly demonstrated by high resolution N-body simulations, 
e.g. [59]. Most often this is taken into account via an effective overall boost factor that 
multiplies the fluxes (see [60] for the earliest analyses). In reality, however, the phenomeno- 
logical implications of substructure are more complicated than that. Indeed, the intensity 
and morphology of the DM annihilation signal is highly sensitive to the way the substruc- 
ture mass function and the subhalo concentration parameters are extrapolated down to 
several orders of magnitude below the actual resolution of the numerical simulations. As 
an effect of propagation and, as a consequence, of the different galactic volumes that con- 
tribute to the signal at Earth, the boost factor can be energy dependent and is in general 
different for different species (e.g. positrons vs antiprotons) [61]. Moreover, in models in 
which the DM annihilation rate is enhanced by the Sommerfeld effect (see [62, 63, 64], and 
then [65, 66]) there are claims that the contribution from substructures might outweigh 
the DM annihilation signal from the smooth main halo alone, e.g. [67]. The gamma ray 
signal of DM annihilation from the Galaxy with a simply modeled population of subhalos 
is presented in [68]. This is however still an active field of research and before the situation 
is somewhat better established we feel it too early to complicate the calculations with a 
highly uncertain additional component. 



3 Fluxes at production 

We consider DM annihilations (parameterized by the DM DM cross section av) and decays 
(described by the DM decay rate Y — 1/r) into the following primary channels: 

e L e L' e %, e Ri ^L^L' ^R^R' T L T Li T R T Ri 

qq, cc, bb, it, 77, gg, 

W+W£, W+Wt, Z l Z l , Z t Z Tj 

hh, 

V e V e , V T V T , 

VV^Ae, W->4fi, VV ->■ 4r, 

where q = u,d,s denotes a light quark and h is the Standard Model (SM) Higgs boson, 
with its mass fixed at 125 GeV [69]. The last three channels denote models in which the 
annihilation or decay first happens into some new (light) boson V which then decays into 
a pair of leptons, along the lines of the models proposed in [65, 70]. 
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The particles produced in Dark Matter annihilation/decay will be provided with parton 

showers and hadronization, in such a way to obtain the fluxes of e^,p, d, 7, ^ e ,n,T at the 
production point. To this goal, we shall use the two most widely used Monte Carlo simula- 
tion programs: Pythia [71] (version 8.135), already used in most DM studies carried out 
so far, and Herwig [72] (version 6.510). In fact, the algorithms implemented in Herwig 
and Pythia are quite different, in both parton showers and hadronization, which makes 
compelling the employment of both codes for the sake of comparing and estimating the 
Monte Carlo uncertainty on a prediction. 

Polarizations and EW corrections. Before moving forward, a brief discussion on po- 
larizations and ElectroWeak corrections is in order. In their current versions, Herwig 
and Pythia contain lepton- or P^-pair production processes in resonance decays, but, un- 
like the channels listed in eq. (2), leptons and vector bosons are treated as unpolarized. 
Furthermore, in the terms which will be clarified later, parton showers include gluon and 
photon radiation, but not the emissions of W's and Z's. However, as discussed in [73] (see 
also [74]) electroweak radiation effects can be particularly relevant for the leptonic and 77 
channels. In fact, the emission of W's and Z's leads to further hadrons in the final state, 
and therefore it significantly modifies the flux of 7's and at energies E <C M, M being 
the DM mass. Moreover, W/Z radiation leads to a p contribution, which is instead absent 
if weak corrections are neglected; this is also true for the the neutrino channels, that also 
yield e s, 7's and p's. In fact, W/Z radiation results in contributions enhanced by one or 
more powers of ln(M/M\yz), with M 3> Mw t z, which do not depend on the DM model. 
Also, in the EW radiation processes, the polarization of the leptons (Left- or .Right-handed 
fermion) and of the massive vectors (Transverse or Longitudinal) plays a role. 

Electroweak emissions can be added to the Pythia event generator following the lines 
of [73], wherein one accounts for the logarithmically-enhanced contributions due to W/Z 
radiation, at leading order in the electroweak coupling constant, as well as leptons/vector- 
bosons polarizations. Of course, the corresponding unpolarized channels can be recovered 
by means of the following averages: 

e+e - = gK + ^gj , w+w . = 2W+Wt+W+W^ 
2 3 

and analogously for r + r~ and ZZ. As electroweak radiation has not been imple- 

mented yet in Herwig, when comparing the two Monte Carlo codes and presenting the 
primary fluxes in section 3.2, we will assume unpolarized particles and turn W/Z emissions 
off. Later on, we will employ only the Pythia code for an extensive study of the par- 
ticle fluxes, and therein the contribution of the large leading electroweak logarithms will 
be taken into account. These effects of W/Z emissions are also included in the numerical 
results collected in the website [29]. 

Let us now go back to the discussion of the list of primary channels in eq. (2). Our 
approach is to consider all these channels on an equal footing, in a manner which is com- 
pletely independent of the DM model. In any given model, annihilation or decay branching 
ratios into the specific channels will instead be dictated by the underlying theory. Some 
channels (such as 77, uu, gg) are 'unusual' as they are often suppressed in many models, 
but from a model-independent point of view they are as viable as any other, so that we 
shall include them and discuss them further below. Operationally, as discussed e.g. in [64], 
s-wave non-relativistic DM DM annihilation can be seen as equivalent to the decay of a 
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Qi resonance with mass M a = 2Mdm, where Mdm is the DM particle mass. Decays of 3i 
into any pair of Standard Model particles can therefore be computed and implemented in 
Monte Carlo generators, along the lines which will be described hereafter. 

The annihilation into SM Higgs, tau, photon and gluon pairs deserves further comments. 
hh. Now that a particle consistent with the SM Higgs boson has been discovered [69], 
we need more than ever to include the corresponding channel in the list of possible an- 
nihilations, and indeed we do so. On the other hand, the detailed properties of such a 
particle are obviously still under very active investigation, so that we have to make some 
guesses/assumptions. For its mass, we assume m h = 125 GeV. For its branching ratios, 
we take those predicted by the Standard Model and embedded in the MonteCarlo codes. 
The values in Herwig and Pythia can differ by up to 25% for a light higgs: Herwig 
has a slightly smaller BR into WW and ZZ with respect to Pythia, while it has has a 
slightly larger BR into bb. Such discrepancies are due to the different accuracy which is used 
to compute the partial widths (see, e.g., [75]). For example, in the decays into WW/ZZ 
Pythia allows both vector bosons to be off-shell, whereas in Herwig at least one is forced 
to be on-shell. In the rate of h —> bb processes, Herwig includes also the resummation of 
mass logarithms ~ ag(m^) ln n (m/ l /m&), which are not resummed in Pythia. Hereafter, we 
shall stick to the default branching fractions for the two codes: more accurate results can 
be obtained by means of the Hdecay program [76], whose use is nevertheless beyond the 
scopes of the present paper. Above all, we stress that the mentioned branching ratios are 
obtained for the Standard Model Higgs boson and that, Beyond the Standard Model, the 
Higgs decay fractions will clearly be different. Should the investigations at the LHC high- 
light a non-SM behavior of the h particle, these assumptions will clearly have to be revised. 
For BSM scenarios both Herwig and Pythia are able to read external data files, provided 
e.g. by the Isajet/SuGra [77] package, containing masses, lifetimes and branching ratios 
of the particles predicted by the chosen new physics model, including non standard Higgs 
bosons. A study of Dark Matter annihilation and decay in scenarios such as the MSSM is 
therefore feasible along the same lines as our current study, but in this work we prefer not 
to make any hypotheses on possible Dark Matter models and we shall stick to decays only 
into Standard Model particles. 

t + t~. As for r leptons, Herwig and Pythia treat them as unpolarized and implement 
the Standard Model three-body decay matrix elements. Alternatively, the two Monte Carlo 
codes can be interfaced with the Tauola package [78], which fully includes polarization ef- 
fects and implements several lepton and hadron decay modes, by means of hadronic matrix 
elements. In the following, we shall nonetheless use the standard Herwig and Pythia rou- 
tines even for the purpose of r decays and subsequent showers and hadronization. In fact, 
this is a reasonable approximation for the observables which we shall investigate, namely 
the hadron/lepton/photon energy fraction in the Dark Matter rest frame and averaged over 
many, many events. A remarkable impact of the inclusion of the r polarization should in- 
stead be expected if one looked at other quantities, such as angular correlations between the 
t decay products from the same event. Furthermore, for the sake of consistency, we prefer 
to use everywhere in our analysis the modelling of hadronization contained in Herwig and 
Pythia, rather than the non-perturbative matrix elements incorporated in the Tauola 
package. 

77. We include 77 as a primary channel: Dark Matter, being dark, has no tree-level cou- 
pling to photons, but 77 production can occur at one loop. This is not to be confused 
with photons emitted by charged particles or produced in three-body annihilations or ra- 
diative hadron decays, such as tt° — > 77. Photons in final-state showers or hadron decays 
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are of course included in the fluxes yielded by Herwig and Pythia (see below for more 
details). Including instead DM annihilation into three-body final states would require a 
specific model of Dark Matter (see e.g. [79, 80]), whereas in this work we shall stick to 
model-independent results. 

gg. Neglecting the case of colored Dark Matter, the DM DM — > gg mode can also take 
place only at one loop. In the Monte Carlo codes which will be employed later on, we shall 
implement the *2s — >■ gg decay in the same fashion as h — > gg, i.e. with an effective @gg 
vertex, assuming that DM is color neutral. 

In the following, we shall first review the basics of parton cascades (section 3.1), then 
we will compare Herwig and Pythia fluxes for a few channels and fixed Mdm (3.2) and 
finally present the Pythia fluxes for several modes and values of the DM mass (3.3). 



3.1 Parton shower algorithms 

In this subsection we review the basics of the parton shower algorithms in Herwig and 
Pythia, focussing in particular on the differences, with the aim of gaining some insight in 
the interpretation of possible discrepancies in the predictions of the two codes. 

QCD (quark, gluon) final state radiation: Showering algorithms rely on the uni- 
versality of the elementary branching probability in the soft or collinear approximation. 
Referring first to quark/gluon final state radiation, the probability to radiate a soft or 
collinear parton reads: 

In (3), P{z) is the Altarelli-Parisi splitting function, z is the energy fraction of the radiated 
parton with respect to the emitter, Q 2 is the evolution variable of the shower. In Herwig 
[72], Q 2 is an energy- weighted angle 4 , which corresponds to angular ordering in the soft 
limit [81]. In the traditional Pythia [82] algorithm, implemented in fortran language, Q 2 
is the virtuality of the radiating parton, with an option to veto branchings that do not fulfil 
the angular ordering prescription. Moreover, the latest version offers, as an alternative, 
the possibility to order final-state showers according to the transverse momentum (kr) of 
the emitted parton with respect to the parent one [83]. kx is also the evolution variable 
in the novel object-oriented Pythia 8 generator [71], written in C++. In either options, 
the Pythia evolution variable is not completely equivalent to angular ordering in the soft 
limit: although in several cases the actual ordering does not make really big changes, 
when comparing with experimental observables sensitive to colour coherence, as done e.g. 
in Ref. [84], Herwig agrees with the data better than Pythia. fc^-ordering yields a 
better inclusion of angular ordering, but nevertheless, as discussed in [85], discrepancies 
with respect to Herwig are still present when considering, e.g., the so-called non-global 
observables [86], sensitive to the radiation in a limited part of the phase space, such as the 
transverse energy flow in a rapidity gap. Throughout this paper, we shall use the Pythia 
8 shower modelling, i.e. ordering in transverse momentum. 

In (3) /S. S (Q\,Q\) is the Sudakov form factor, expressing the probability of evolution 
from Q\ to Q\ with no resolvable emission. In diagrammatic terms, it sums up all virtual 

4 In the Herwig showering frame, Q 2 ~ E 2 (l — cos 8), where E is the energy of the splitting parton 
and is the emission angle [81]. 
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and unresolved real emissions to all orders. In particular, the ratio of form factors in (3) 
represents the probability that the considered emission is the first, i.e. there is no emission 
between Q 2 and Qm ax , where Q^ax i s se ^ by the hard-scattering process. Ql is instead the 
value of Q 2 at which the shower evolution is terminated and hadronization begins. The 
Sudakov form factor is given by the following equation: 

. f [ Q ™*dk 2 , a s (z,k 2 ) n/ A 

As(QL x , Q ) = exp <^ - / — / dz ^ P(z) \ . (4) 

In Eq. (4), the lower integration limit is related to the evolution variable Q and the shower 
cutoff Qq: it is given by z m[n = Q\jQ 2 in Pythia and z m i n = Qo/Q in Herwig, consistently 
with virtuality /angular ordering. The upper limit is instead z max = 1 — z min . For any Q, 
such conditions imply a larger ^-evolution range in Pythia with respect to Herwig. 

For multiparton radiation, iterating the branching probability (3) is equivalent to per- 
forming the resummation of soft- and collinear-enhanced radiation. As discussed, for ex- 
ample, in [87] in the framework of the Herwig event generator, parton shower algorithms 
resum leading logarithms (LL) in the Sudakov exponent, and include a class of sub leading 
next-to-leading logarithms (NLL) as well. The strong coupling constant in (3) is evalu- 
ated at the transverse momentum (k T ) of the radiated quark/gluon with respect to the 
parent parton: in this way, one includes in the showering algorithm further subleading 
soft/collinear enhanced logarithms [87]. Moreover, in Herwig the two-loop coupling con- 
stant in the physical CMW scheme [87] is used 5 ; in Pythia the /3 function is instead 
included in the one-loop approximation. The actual value of as{Mz) comes after fits to 
LEP data, which yield a s {M z ) ~ 0.116 in Herwig and a s {M z ) ~ 0.127 in Pythia. 
Different values of as[Mz) would lead to different total cross sections for hard-scattering 
processes mediated by the strong interaction, but have very little impact on differential dis- 
tributions, such as the ones which we will investigate, as long as tuned versions of Herwig 
and Pythia are used (see e.g. the discussion in [75]). 

Also for the purpose of hadronization, the two programs implement very different mod- 
els, namely the cluster model [88] (Herwig), based on colour preconfinement and closely 
related to angular ordering, and the string model [89] (Pythia), both depending on a few 
non-perturbative parameters. These parameters, along with other quantities, such as the 
shower cutoff or quark and gluon effective masses, are fitted to experimental data, e.g., 
e + e~ or pp data from LEP and Tevatron experiments (see, e.g., Refs. [90, 91]). In prin- 
ciple, whenever one runs Herwig or Pythia at much higher energies, as will be done 
in the following, such fits may have to be reconsidered. In fact, although the hadroniza- 
tion transition is universal, when using a hadronization model along with a perturbative 
calculation or a parton shower algorithm, one assumes that the hadronization model even 
accounts for the missing perturbative contributions, which are clearly process-dependent. 
Therefore, whenever one has data from other experiments, one should check that the best- 
fit parametrizations are still able to reproduce the data. In this paper, for the sake of a 
consistent comparison, we shall use the default values of the parameters employed in the 
default versions of Herwig and Pythia, which were fitted to LEP data (i.e. at the typical 
100 GeV energy scale). In the future, data at higher energies (e.g. from the LHC) will 

5 The CMW scheme consists in denning a Monte Carlo strong coupling a^ c related to the usual MS 

one via a^ c = ag IS [l + Ka^ s /(2tt)], where the explicit expression for K can be found in [87]. By means 
of this rescaling, the Herwig Sudakov form factor includes threshold-enhanced corrections in the NLL 
approximation. 
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be available and will possibly allow to retune the two Monte Carlo codes to better mimic 
high-mass Dark Matter annihilation/decay. 

QED (photon) final state radiation: Extending the algorithm (3) to include 
photon radiation off quarks and leptons, as well as photon branching into quark or lepton 
pairs, is straightforward. Pythia (the version 8.135 that we use) include all such processes, 
but it does not include photon radiation from W + W~ final states. We have therefore to add 
this by hand as discussed in [73] (see footnote 2 in vl on the arXiv). The latest version of 
HERWIG, on the other hand, contains q — > branchings, but does not implement photon 
emissions off leptons as well as 7 — > f f processes, which are instead present in Herwig++ 
[92], the object-oriented version written in C++, along the lines of [93]. In the following 
we will nevertheless employ the latest fortran version of Herwig and delay to future work 
the use of HERWIG++. We can indeed anticipate that, in the phase space regions in which 
we will mostly be interested, the partial lack of QED radiation in fortran Herwig has a 
small effect in the comparison with Pythia. Also, we should always keep in mind that, 
although shower and hadronization are implemented in a different fashion, the two codes 
have been fitted to agree with LEP data (as discussed above), and therefore possible lacks 
are somehow compensated by means of suitable parameter tuning. 6 



3.2 Dark Matter fluxes comparing Herwig and Pythia 

We present the energy spectra of final-state particles from DM annihilation, yielded by 
Herwig and Pythia. For this purpose, as anticipated, we modified the two codes in order 
to include the decay of a generic resonance — > ab, in such a way that we are allowed to 
fix the mass M @ = 2M and specify particles a and b. Moreover, we changed the hadron 
decay tables, to make it sure that hadrons, such as kaons or pions, which are often treated 
as stable when simulating collider phenomenology, decay according to the branching ratios 
quoted in the PDG [51]. Our results will be expressed in terms of the energy fraction 

K 

where K is the kinetic energy of the final-state stable hadrons/leptons/photons in the rest 
frame of S> '. We shall plot the particle multiplicity as a function of the logarithmic energy 
fraction, i.e. dN/dlogx; our spectra will be normalized to the average multiplicity in the 
simulated high-statistics event sample. Also, as pointed out before, this comparison will 
be carried out for production of unpolarized particles and without including any effect of 
final-state weak boson radiation. 

An example of the comparison of the DM fluxes from Pythia and Herwig is presented 
in Fig. 2, where we show the photon, electron, antiproton and neutrino dN/dlogx spectra 
for the channels DM DM — > qq, gg, W + W~ and t + t~. In Fig. 2 we have set the DM mass 
to Mdm = 1 TeV, but we can anticipate that similar dN / d log x hold for all DM masses 
.Mdm 3> Mz,m t . Astrophysical experiments are currently probing i\~<100GeV, whose 
corresponding range of x depends on the chosen Mdm; in particular, the low-x tails mostly 
determine the DM signals if Mdm is very large. Overall, we note the following features: 

• For the qq modes there is a reasonable agreement between Pythia and Herwig, 
for all final-state particles and through the whole x spectrum, including the low- 



6 We also point out that, although in future perspectives the C++ generators will supersede the fortran 
ones, fortran Herwig and Pythia will still be used as they provide parton shower and hadronization to 
the so-called 'matrix-element generators' like Alpgen [94] or Madgraph [95]. 
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Figure 2: Comparison between Monte Carlo results: Pythia is the continuous line, Her- 
WIG is dashed. Photons (red), e^ (green), p (blue), v = v e + + v T (black). 



energy tails. In fact, although the centre-of-mass energy has been increased to 2 
TeV, the S — y qq is similar to Z/j* — y qq processes at LEP, which were used when 
tuning the Herwig and Pythia user-defined parameters. Nevertheless, we note some 
discrepancy, about 20%, especially in the neutrino spectra, as Pythia yields overall 
a higher multiplicity, and in the p distribution, where Herwig is above Pythia 
especially at large x. 

• Some discrepancy, up to a factor of 2, is instead found for the gg mode (which is, 
however, presumably not the dominant one in DM phenomenology). In fact, unlike 
the qq mode, the S —y gg channel does not have a counterpart at LEP; the differences 
in parton showers and hadronization in Herwig and Pythia, as well as the fact that 
we are running the two codes at a much higher energy with respect to LEP, may thus 
be responsible for this discrepancy. In detail, as far as the 7, e ± and p spectra are 
concerned, Herwig is above Pythia at small x and below at large x; the Pythia 
neutrino multiplicity is instead above the Herwig one in the whole x range, especially 
for x > 1CT 5 . 

• Lepton modes (here exemplified by the t~t + case) exhibit a significant disagreement, 
especially in the photon spectra, where Pythia yields a remarkably higher multi- 
plicity with respect to Herwig for x < 10~ 2 . As we pointed out before, Pythia 
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includes £ — > £7, £ being a charged lepton, and 7 — > ff processes, which are instead 
absent in HERWIG, whose photons may come only from q — > or radiative hadron 
decays. The lack of this type of 7 radiation in Herwig might therefore explain the 
discrepancy in the photon spectra in the t + t~ channel. 

• The fourth panel shows an example of 'composite' mode, i.e. the DM DM — > W + W~ 
channel. This process, where a W pair is produced via annihilation of colourless 
particles, is alike e + e~ — > W + W~ at LEP 2, where Herwig and Pythia have been 
carefully tested. In fact, the electron and positron spectra exhibit good agreement, 
and even the antiproton multiplicities are not too far, although Herwig yields a 
higher dN/dlogx for x > 10~ 3 . More visible discrepancies are instead exhibited by 
the neutrino spectra, as the Pythia multiplicity is higher than Herwig throughout 
all x range, and by low-x photons, where Pythia is much above Herwig. The latest 
discrepancy can be traced back to the differences in photon radiation off leptons 
commented on above. 

Similar features are observed for other values of the DM mass and for other modes. 
In view of these considerations, mainly the fact that Pythia includes also £ — > £7 and 
7 — > ff processes and its spectra can be supplemented by large electroweak logarithms 
according to [73], hereafter we shall employ only Pythia to present the spectra for several 
decay/annihilation modes and mass values. The Monte Carlo uncertainty, gauged from 
the comparison between Pythia and Herwig, can be estimated to be on average about 
±20%, with the exception of the (probably negligible) gluon-gluon channel and photon 
(neutrino) spectra in the lepton (VTV4 7 ) mode, which instead exhibit a larger disagreement, 
especially for very small values of x. As pointed out above, this latest discrepancy will 
likely get milder if one implemented further QED-type branchings in Herwig or possibly 
used HERWIG++. 

3.3 Dark Matter fluxes: results 

We therefore compute the fluxes ' of e^, p, d, 7, v e /vr in a large range of DM masses Mdm = 
5 GeV — > 100 TeV, by using the Pythia event generator, and provide them in numerical 
form on the website [29], both in the form of Mathematica® interpolating functions and 
numerical tables. 8 Such computing-power demanding results have been obtained using the 
EU Baltic Grid facilities [96]. 

In fig. 3 we present some examples of the spectra produced by the annihilation of two 
DM particles with mass Mdm (normalized per annihilation), for four values of Mdm- They 
correspond to the fluxes from the decay of a DM particle with mass 2Mdm- 

Some specifications on these fluxes are in order. About all fluxes: The fluxes presented 
from here on and employed in the rest of the paper include EW corrections, as discussed 
above and in [73]. However, in [29] we also provide, for comparison, all the spectra before 
EW corrections. 

About 7 ray fluxes: We specify that of course the fluxes here include only the prompt 
emission and not the secondary radiation (e.g. due to Inverse Compton processes) that 
we discuss in Sec. 6. Furthermore, we recall that by 'prompt emission' we here mean all 

7 We note that a similar work has been performed recently, limited to gamma rays, in [97]. 
8 Analogous results yielded by the Herwig code can be obtained by contacting the authors (G.C. and 
F.S.). 
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Figure 3: Primary fluxes of e , p, d, 7 and u e . 
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Figure 4: Energy distribution between the final states particles: e , hadrons (p + d), 
7 and for a set of characteristic annihilation channels. The inner ( outer) pie refers to a DM 
mass of 200 GeV (5 TeV). For each pie chart, the first caption gives the energy fraction going 
into 7 and e^ (E 1+e ) with respect to the total. The second caption gives the energy fraction into 
hadronic final states (E p+( i) with respect to 7 and e^. 
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photons in final-state showers or hadron decays as given by Pythia, including those from 
(IR-enhanced model-independent) QED and EW bremsstrahlung as discussed above and 
in [73]. But further contributions to prompt emission can come from other three-body 
final states such as internal bremsstrahlung [79, 80]: these can only be computed in the 
framework of a precise DM model because one needs to know the higher order QED anni- 
hilation/decay diagram. These are not included. 

About fluxes of anti-deuterons: They are computed taking into account the jet structure of 
the annihilation products scale with the cube of the uncertain coalescence parameter, here 
fixed to po = 160 MeV; for details on the computation we refer the reader to [98] . 
About fluxes of neutrinos and anti-neutrinos: Those that we provide here are of course the 
neutrino spectra at production; the corresponding fluxes at detection are affected by oscil- 
lations (if travelling in vacuum, such as for neutrinos from DM annihilations/decays at the 
Galactic Center) and/or by interactions with matter (if e.g. from DM annihilations/decays 
in the center of the Sun). The fluxes at detection of neutrinos having traveled in vacuum 
from a distant astrophysical source can be obtained taking into account average oscillations 
with the formula 

3 / 0.6 0.2 0.2 \ 

P{y t v v ) = P{v t v t ) = ~ °- 2 °- 4 °- 4 (6) 

i=i \ 0.2 0.4 0.4 J 

where i runs over neutrino mass eigenstates, and the elements \V^\ of the neutrino mixing 
matrix can depend on its unknown CP-violating phases. The case of neutrinos from the 
center of the Sun is more complicated: they can be obtained along the lines of the analysis 
in [99] and we leave to future work their detailed presentation. Finally, we recall that 
neutrinos detected after having crossed the Earth can experience additional oscillation 
effects. 

To conclude this quantitative presentation of the DM fluxes, in fig. 4 we show some 
characteristic energy distributions between the final-states particles: e , hadrons (p + d), 
7 and v. The inner pie refers to a DM mass of 200 GeV (a typical SuSy WIMP value) 
while the outer pie to 5 TeV (taken as a typical multi-TeV case). One can see that the 
portion of energy which goes into gamma rays and e is often the most important one and 
always dominates over the energy fraction of the hadronic final states for all the channels. 
This is especially relevant in the context of extragalactic gamma rays signatures, where the 
energy fraction in is quickly converted to gamma rays due to Inverse Compton radiation. 
For the channels involving and t + t~ and of course for the vv channels, the portion 

of energy carried away by neutrinos becomes the dominant one. The fractions are rather 
independent of the mass of the DM particle, with some exceptions. For example, in the 
vv channels, primary neutrinos start to radiate gamma rays and charged leptons due to 
radiative weak corrections when Mom is above the electroweak scale (i.e. for the outer pie 
in the figure) and this increases the energy fraction of 7 and e . 

4 Charged Cosmic Rays 

Having at disposal the energy spectra of charged particles per annihilation at production, 
as generated by MonteCarlos and as discussed in the previous section, we next need to 
consider where these fluxes of particles are produced in the galactic halo and how they 
propagate to the Earth. 
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For simplicity, we present separately the propagation formalism for electrons or positrons, 
for antiprotons and for antideuterons. In this latter case, only a few trivial changes have 
to be implemented with respect to antiprotons, as we discuss below. 

In general one ends up with a convenient form for the propagated fluxes in terms of a 
convolution of the spectra at production with a propagation function that encodes all the 
intervening astrophysics: see eq. (13) for e (or eq. (22) for the approximated treatment), 
eq. (27) for p and eq. (33) for d. 

4.1 Propagation functions 

4.1.1 Electrons or positrons: full formalism 

The differential e flux 9 per unit of energy from DM annihilations or decays in any point in 
space x and time t is given by rf$ e ± / dE (t, x, E) = v e ± f /Arc (units 1 / GeV ■ cm 2 • s ■ sr) where 
v e ± is the velocity (essentially equal to c in the regimes of our interest). The e ± number 
density per unit energy, f(t,x,E) = dN e ±/dE, obeys the diffusion-loss equation [100]: 



with diffusion coefficient function JC(E, x) and energy loss coefficient function b(E, x). They 
respectively describe transport through the turbulent magnetic fields and energy loss due 
to several processes, such as synchrotron radiation and Inverse Compton scattering (ICS) 
on CMB photons and on infrared or optical galactic starlight, as we discuss in more detail 
below. Notice that other terms would be present in a fully general diffusion-loss equation 
for Cosmic Rays, such as diffusive re-acceleration terms (describing the diffusion of CR 
particles in momentum space, due to their interactions on scattering centers that move 
in the Galaxy with an (Alfven) velocity V a ) and convective terms. These are however 
negligible for e ± , see e.g. [100, 101]. Eq. (7) is solved in a diffusive region with the shape 
of a solid flat cylinder that sandwiches the galactic plane, with height 2L in the z direction 
and radius R = 20 kpc in the r direction [102]. The location of the solar system corresponds 
to x = (r ,£ ) = (8.33 kpc, 0). Boundary conditions are imposed such that the e ± density 
/ vanishes on the surface of the cylinder, outside of which electrons and positrons freely 
propagate and escape. 10 Assuming that steady state conditions hold (as it is if one assumes 
that the typical time scales of the DM galactic collapse and of the variation of propagation 
conditions are much longer than the time scale of propagation itself, of the order of 1 Myr 
at 100 GeV energies [104]), the first term of eq. (7) vanishes and the dependence on time 
disappears. 11 Before illustrating the solution method, we briefly comment on the different 
pieces of the equation. 

9 Notice that with the notation e ± we always refer to the independent fluxes of electrons eT or positrons 
e + , which share the same formalism, and not to their sum (for which we use the notation e + + e~ when 
needed) which of course differs by a trivial factor 2. 

10 See [103] for the impact of not neglecting the propagation outside the cylinder. 

U A caveat on this point is that the time-independence of the diffusion process might not be justified in 
extreme environments such as the galactic central regions, where the propagation conditions may possibly 
change on a short enough time scale that they make this assumption invalid. E.g. recently the Fermi 
satellite has pointed out the existence of large gamma-ray structures (dubbed 'Fermi bubbles') above and 
below the Galactic Center [105]. A detailed modeling of the impact of these possible features on CR 
propagation is, for the moment, well beyond the scope of our analysis and probably of most DM related 
ones. 
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Figure 5: Energy loss coefficient function for electrons and positrons in the Milky Way. 
Left panel: at several locations along the galactic radial coordinate r, right panel: above (or below) 
the location of the Earth along the coordinate z. The dot points at the value of t q (see next 
subsection) . 



The e ± energy loss coefficient function b(E, x) is in general dependent on the position 
x, since the energy losses suffered by the e ± are sensitive to the environment: 

- d ^^b(E,x) = ^E 2 u , u = u B (x) + J2^(x)R? N (E), (8) 

e i 

where o~t = 8nr 2 /3, with r e = a em /m e , is the Thompson cross section. The first addend 
in u is associated with synchrotron losses, the second one with ICS losses 12 . ub = B 2 /2 is 
the energy density in galactic magnetic fields B and w 7)i = f dEni(E) is the energy density 
in light. Here i runs over the three main components: CMB, star-light and dust-diffused 
InfraRed light. For the CMB, n(E) is just the black body spectrum with T = 2.725 K 
and one gets m 7i cmb = 0.260 eV/cm 3 . For IR light and starlight, we extract the maps of 
their distribution and energy profile in the Galaxy from Galprop [30] (an approximated 
formalism which employs a superposition of black body spectra also for IR and starlight 

12 So one can also define 

b syn (E,x) = ^E 2 u B (x) and in c (E,g) = ^ E 2 u yii (x)Rf N (E) (9) 
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has been discussed in [106]). The x dependence in b is due to the fact that the composition 
of the background light for ICS is different in different points of the halo (e.g. the center or 
the periphery of the Galaxy) and the value of the magnetic field also varies (much higher 
in the center than elsewhere). 

The dependence of b on the energy E, in turn, is dictated by the energy dependence of 
the rates of the different loss processes. In particular, for IC scattering one has b oc E 2 as 
long as the scattering happens in the Thomson regime, where the factor R™(E) = l. 13 For 
large enough electron energy the IC scattering enters into the full Klein-Nishina regime, 
where the 7e ± scattering rate becomes rapidly smaller than the Thomson approximation, 
and thus Rf N (E e ) < 1 (see e.g. fig. 2 of [107]) reducing b(E). The synchrotron loss rate, 
instead, is always proportional to the square of the electron/positron energy E 2 . 

The profile of the magnetic field in the Galaxy is very uncertain and we adopt the 
conventional one 

B(r, z) = B exp[-(r - r e )/r B - \z\/z B ] (10) 

as given in [108], with B = 4.78 /zG, r B = 10 kpc and z B = 2 kpc. With these choices, 
the dominant energy losses are due to ICS everywhere, except in the region of the Galactic 
Center and for high e energies, in which case synchrotron losses dominate. All in all, 
the b(E,x) function that we obtain is sampled in fig. 5 and given in numerical form on 
the website [29]. In the figure, one sees the E 2 behaviour at low energies changing into a 
softer dependence as the energy increases (the transition happens earlier at the GC, where 
starlight is more abundant, and later at the periphery of the Galaxy, where CMB is the 
dominant background). At the GC, it eventually re-settles onto a E 2 slope at very high 
energies, where synchrotron losses dominate. 

The diffusion coefficient function /C is also in principle dependent on the position, since 
the distribution of the diffusive inhomogeneities of the magnetic field changes throughout 
the galactic halo. However, a detailed mapping of such variations is prohibitive: e.g. they 
would have different features inside/outside the galactic arms as well as inside/outside the 
galactic disk, so that they would depend very much on poorly known local galactic geogra- 
phy. Moreover, including a spatial dependence in K, would make the semi-analytic method 
described below much more difficult to implement numerically. We therefore leave these 
possible refinements for future work 14 and, as customary, we adopt the parameterization 
1C(E, x ) = K, {E/ GeV) 5 = Kq e 5 . 

The values of the propagation parameters 5, K and L (the height of the diffusion 
cylinder defined above) are deduced from a variety of cosmic ray data and modelizations. 
It is customary to adopt the sets presented in Table 1, which are found to minimize or 
maximize the final fluxes. 15 

Finally, DM DM annihilations or DM decays in each point of the halo with DM density 

13 We recall that the Thomson regime in electron-photon Compton scattering is identified by the condition 
e max = < m e , where e denotes the energy of the impinging photon, e' the same quantity in the rest 
frame of the electron, 7 is the Lorentz factor of the electron and m e is the electron mass. When e scatter 
on CMB photons (e ~ 2 10~ 4 eV) the condition is satisfied up to ~ TeV e ± energies. For scatterings on 
more energetic starlight (e ps 0.3 eV), the condition breaks down already above ps few GeV e 4 " energies. 

14 See [109] for a recent analysis for antiprotons. 

15 We stress, however, that the determination of these parameters is a whole evolving research area, which 
will certainly update these values in the future as more refined modelizations and further CR data become 
available. See e.g. [112, 113, 114] for recent references. The choices presented in Table 1 should be seen as 
the current bracketing of sensible possibilities. 
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Model 


J-JlCLy \jL Ullo \JL JJColLIvJllo 

5 /C [kpc 2 /Myr] 


6 /C [kpc 2 /Myr] Konv [km/s] 


L [kpc] 


MIN 


0.55 0.00595 


0.85 0.0016 13.5 


1 


MED 


0.70 0.0112 


0.70 0.0112 12 


4 


MAX 


0.46 0.0765 


0.46 0.0765 5 


15 



Table 1: Propagation parameters for charged particles in the Galaxy (from [110, 111]). 



p(x) provide the source term Q of eq. (7), which reads 

Q= l{jt^) 2f ^' ^^B'*^ (annihilation), (11) 

Q= {iL)fr ^ = ^ T nk (decay) - (12) 

where / runs over all the channels with e in the final state, with the respective thermal 
averaged cross sections av or decay rate T. 



4.1.2 Electrons or positrons: result 

The differential flux ofe ± d§ e ±/dE = v e ±f /Air in each given point of our Galaxy for any 
injection spectrum can be written as 



f 1 /p(f)\ 2 [ m ™ dN f e ._ 



dE ^' x) - Airb(E,x) 



2 V M nM ) y 



/■«DM dN + 

(crv)f / dE s i, (E s ) I(E,E B ,x) (annihilation) 
Je dE 



1 M J2 T f / dE s —f(E s )I(E,E s ,x) (decay) 

f *J E 

(13) 

where E s is the e energy at production ('s' stands for 'source') and the generalized halo 
functions I(E,E s ,x) are essentially the Green functions from a source with fixed energy 
i^s to any energy E. In other words, the halo functions I encapsulate all the astrophysics 
(there is a halo function I for each choice of DM distribution profile and choice of e 
propagation parameters) and are independent of the particle physics model: convoluted 
with the injection spectra, they give the final spectrum searched for. They obey I(E, E, x) = 
1 and I(E,E s ,x) = on the boundary of the diffusion cylinder. Neglecting diffusion (i.e. 
setting /C = 0) one would have I(E,E s ,x) = 1. These functions are provided numerically 
on the website [29] in the form of Mathematica® interpolating functions. Plugged in 
eq. (13), they allow to compute the e flux everywhere in the Galaxy. 

The functions particularized to the location of the Earth, that is: I(E, E e , r ), are 
plotted in fig. 6 and provided numerically on the website [29] too. Plugged in eq. (13), 
these allow to compute the e ± flux at the location of the Earth, $(e, r , z ). We also 
provide separately the resulting fluxes (see the next subsection 4.2.1). 

The generalized halo functions I are computed as follows (the uninterested reader can 
skip the rest of this section). Due to numerical issues it is convenient to search for the 
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Figure 6: Generalized halo functions for electrons or positrons, for several different val- 
ues of the injection energy eg (color coded). The superimposed dotted black lines are the 'reduced' 
halo functions discussed in 4-1-3. 



solution of eq. (7) using an ansatz similar to eq. (13) but somewhat different: 

f 1 / p ^ 2 <- M ™ 



6 T (e) 



2 VM DM 

Pq 



de B f^(e s )I(e,e s ,x) 



Mi 



DM 



Mdm/2 



de s /^ c (e s )/(e,e s ,f) 



(annihilation) 



(decay) 



(14) 



where e = E/GeV. Here we adopted the (arbitrary but convenient) normalizing factor 
br(e) = e 2 GeV/r , with r = GeV/6(l GeV, x G ) = 5.7 x 10 15 sec, which is the energy 
loss coefficient at Earth in the Thomson limit regime. Plugging now the ansatz (14) in the 
differential equation (7) one can recast (7) into a partial differential equation for I(e,e s ,x) 
(this extends the solution method first discussed (to our knowledge) in [115]). Indeed, (7) 
becomes 



M DM (M DM /2) 



O'0< 
M DM (M DM /2) 



de s f ini {e s ) V 2 /(e,e s ,f) + 



b{e,x) f ini (e)I(e,e s ,x) 



de B f(e 



d fb(e,x) ~ 



at \ b T (e) 



/inj(e) 



(15) 



where r\ = 1, 2 for decay or annihilation scenarios respectively and the upper integration 
limit changes accordingly. One then extracts the partial differential equation for /: 



V 2 /(e,e s ,x) + 



1 d f b(e,x) - 
/C o r e^9 e ^ T (6) ne ' es ' X ' 



0. 



(16) 



with boundary conditions 



&T(e S ) / 



-^(^) ^sj ■'■max) 



6(e s ,f) V Pe J (17) 
= 0, with x max = {R,L). 
Finally the halo functions with the normalization conventions of eq. (13) are obtained as 

(18) 

Solving numerically eq. (16) with (17) allows to compute the I(e,e s ,x) and in turn the 
I(e,e s ,x). 



I(E,E s ,x) = I(e,e s ,x) 





(f)l 


b(e, x) 





4.1.3 Electrons or positrons: approximated energy loss 

The above treatment is pretty general in that it allows to compute the propagated fluxes 
taking into account the full energy and position dependance of b(E, x), as discussed above. 
An approximated formalism had been adopted in the past (see e.g. [110] and references 
therein) and we report it here for completeness and to compare with our full result. 

Assuming a space-independent b = b T (e) = e 2 GeV/r everywhere in the Galaxy, one 
can define a 'reduced' halo function T(\d,x) (and a simplified differential equation for 
it) in terms of a single quantity \d = Ai)(e,e s ) = ^4/Cqt (e 5-1 — ef -1 ) /(l — 5), which 
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Figure 7: Reduced halo function X(A/j) for from 
annihilating Dark Matter, for the different DM profiles 
and sets of propagation parameters, and the corresponding 
fit parameters to be used in eq. (23) . 
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Figure 8: Reduced halo function I(Xd) for e^ from 
decaying DM, for the different halo profiles and sets of 
propagation parameters, and the corresponding fit parame- 
ters to be used in eq. (23). 
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represents the diffusion length of e ± injected with energy es and detected with energy e. 
One has (instead of the full equation (16)) 



V 2 Z(A D ,f) - ^--^-Z(A D ,f) =0, 
Ad a ad 




p{x) 
Pq 



(19) 



Solving eq. (19) provides T{\d,x) in any given point x. 16 With that, one can write the 



16 Alternatively, one can find the solution for X(Xd, x) via an expansion in Bessel and Fourier series [110]. 
Formally one finds that 



T(X D ,r,z)= ^2 Ja(( n r/R) sm (^"( z + L )) ex P 




(20) 



where Ji is the Bessel function of the first kind (cylindrical harmonic) of order i, £ n is the n-th zero of the 
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equivalent of eq. (13) as 



d$ e ± , 1 v e ± 

(e,x) 



dE 



47r & T (e) 



Pq 

■2 V M DM 

Po 



A/dm 



^e s ) I( Ad (e, e s ) , x) (annihilation) 



Mi 



DM 



/;■ / de s —jj^(e s )l(\ D (e,e s ),x) 



(decay) 
(22) 

The function /(Ad, ">"©, £©) at the location of the Earth is well reproduced in terms of the 
fit 



Z{^d) = a + ai tanh 



bx-i 



r/ocx]) ( ] +a 3 



(23) 



with £ = log 10 (A/)/kpc) and the coefficients given in the tables in Figure 7 and 8 and also 
reported on the website [29]. 17 

These approximated halo functions are also superimposed (black dotted lines) in the 
plots of fig. 6, having chosen €5 = 1 GeV and the annihilation case. It is evident that 
the reduced halo functions miss much of the richer structure of the full ones, especially for 
MED and MAX parameters and for peaked profiles. On the other hand, the zeroing of the 
functions in the two approaches occurs at similar values of x in each plot (compare the 
black dotted lines and the orange ones corresponding to es — 1 GeV in the color coding). 
In computing the convolution as an integral over €5, the different shapes and the similar 
zeroing somewhat counteract one each other, so that the final spectra as computed by the 
full formalism will be different, but not drastically, with respect to the approximated result. 
We illustrate later in specific examples (see fig. 13) the quantitative impact. In points of 
the Galaxy other than the location of the Earth, e.g. close to the Galactic Center, the 
difference between the two computations would however be more important. This affects 
the spectra of ICS gamma rays produced from these regions. 



4.1.4 Antiprotons 

The propagation of antiprotons through the galaxy is described by a diffusion equation 
analogous to the one for positrons. Again, the number density of antiprotons per unit 
energy f(t, x, K) = dNp/dK vanishes on the surface of the cylinder at z — ±L and r = R. 
K = E — irip is the p kinetic energy, conveniently used instead of the total energy E (a 

i = function and R n , m corresponds to the Bessel- and Fourier-transform of (p/po) 2 '- 

Rn,m= Jl(c 2 )2 E 2 J drrJ (C n r/R)^J + ^dzahx(mn(z + L)/2L^(^^^ . (21) 

This method, however, converges much more slowly than numerically solving eq. (19), and it is therefore 
unpractical for computing the solution in more than one point (r, z) . We checked anyway that the two 
solutions coincide. In the numerical computations performed with this Bessel-Fourier expansion method it 
is convenient to smooth out the steepness of some of the profiles close to the GC, adopting the prescription 
discussed in [116]. It simply amounts to replacing the divergent profile by a well behaved one below an 
arbitrarily chosen critical radius of r crit = 0.5 kpc from the GC, while approximately preserving the absolute 
number of annihilations in that region. Such well behaved profiles are plotted in fig. 1 as dashed lines. 
We just cut at r cr it/10 for the Einasto profiles. It can be checked that such smoothing does not affect the 
determination of I(A£>, f©, Zq). 

17 The fit functions reproduce the results of our numerical calculation to better than 5%, with the ex- 
ception of the EinastoB case, for which the accuracy drops to a still acceptable 10%, over the whole range 
Ajj = 0.1 — >• 100 kpc. The fit functions should not be used outside of this range. 
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distinction which is of course not particularly relevant when one looks at fluxes originating 
from TeV-scale DM, i.e. at energies much larger than the proton mass m p , but important 
for the low energy tails and in the case of small DM masses). Since m p ^> m e one can 
neglect the energy loss term that was instead important for positrons. But new terms 
appear in the diffusion equation for /, which reads 

df d 

-± - K(K) ■ V 2 f + — (sign(z) / Konv) =Q-2h 8{z) (r ann + r non _ ann )/, (24) 



where: 



The pure diffusion term can again be written as )C(K) = /Co/3 (p/ GeV) 5 , where 

p = (K 2 + 2m p K) 1 / 2 and = v p /c = (l — m 2 /(K + rrip) 2 ) 1 ^ 2 are the antiproton 
momentum and velocity. S and /Co are given in Table 1. 

The V COQV term corresponds to a convective wind, assumed to be constant and directed 
outward from the galactic plane, that tends to push away p with energy T < 10 m p . 
Its value is given in Table 1. 

The source term Q due to DM DM annihilations or DM decay has a form fully 
analogous to eq. (11) or (12), with E now formally replaced by K. 

The first part of the last term in eq. (24) describes the annihilations of p on interstellar 
protons in the galactic plane (with a thickness of h — 0.1 kpc <C L) with rate r ann = 
(tih + 4 2//3 nHe)Cp£ n t>p, where nn ~ 1/cm 3 is the hydrogen density, % e ~ 0.07 wh is the 
Helium density (the factor 4 2//3 accounting for the different geometrical cross section 
in an effective way) and a^p n is given by [117, 118] 



aim 



661 (1 + 0.0115 K-°- 7U - 0.984 /C - 0151 ) mbarn, for K < 15.5 GeV 



a PP \ 36 K~ - 5 mbarn, for K > 15.5 GeV 

(25) 

The second part, similarly, describes the interactions on interstellar protons in the 
galactic plane in which the p's do not annihilate but lose a significant fraction of 
their energy. Technically, one should keep them in the flux, with a degraded energy: 
they are referred to as "tertiary antiprotons" . We adopt instead the simplifying 
approximation of treating them as if they were removed from the flux. The cross 
section that we need for the whole last term of eq. (24) is then the sum of a^ n + 

a non-ann = ^inel ft j g giyen j n [ n? ] ag 

^\ K ) = 24.7 (1 + 0.584 R- - 115 + 0.856 R-°- 5m ) mbarn (26) 

(at large energies this expression has to be replaced by a better approximation [119]). 
We find, anyway, that the precise expressions adopted for these cross sections do not 
significantly impact the final results. 

We neglect, as just said, the effect of "tertiary antiprotons". It can be re-included in 
terms of an absorption term proportional to a different u 11011 - 1111 ^ and of a re-injection 
term Q tert proportional to the integrated cross section over f(K). The full solution of 
the resulting integro-differential equation can be found in [119]. The effect of tertiaries 
is mainly relevant at low energies K < few GeV. 
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- Finally notice that we do not include diffusive reacceleration. It does not play a major 
role for p (as it was not for e ) at the high energies in which we are mostly interested, 
say larger than tens of GeV. It can instead affect the spectrum at ~GeV energies, and 
it can be reintroduced in an effective way by adding an effective energy-loss coefficient 
and/or modifying the energy dependance of the spatial diffusion coefficient. We refer 
to [111, 119, 120] for details. 

Assuming again steady state conditions the first term in the diffusion equation vanishes, 
and the equation can be solved analytically [120, 121, 122]. In the "no-tertiaries" approxi- 
mation that we adopt, the solution for the antiproton differential flux at the position of the 
Earth d^p/dK (K,r®) = Vp/(4Tr)f acquires a simple factorized form (see e.g. [Ill]) 



d$p 
dK 



Pe 
Mi 



DM 



, v dN P 



(annihilation) 



An 



teHX>§ (decay) 



(27) 



The / index runs over all the annihilation channels with antiprotons in the final state, with 
the respective cross sections or decay rates; this part contains the particle physics input. 
The function R(K) encodes all the astrophysics of production and propagation. 18 There is 
such a 'propagation function' for annihilations and for decays for any choice of DM galactic 
profile and for any choice of set of propagation parameters among those in Table 1. We 
provide R{K) for all these cases in terms of a fit function 



log 10 [R(K)/Myr] 



do + Ctl K + Ct2 K 2 + (2 3 K 3 + a^K 4 + CI5 ft 5 , 



(30) 



with k = log 10 Kj GeV and the coefficients reported in the tables in fig. 9 and 10 (and also 
reported on the website [29]). 19 

Finally, for completeness we also mention the average solar modulation effect, although 
it is mainly relevant for non-relativistic p's: the solar wind decreases the kinetic energy 
K and momentum p of charged cosmic rays such that the energy spectrum dQp^/dK^ of 
antiprotons that reach the Earth with energy and momentum (sometimes referred to 
as Top of the Atmosphere 'ToA' fluxes) is approximatively related to their energy spectrum 
in the interstellar medium, dQp/dK, as [123] 



dK a 



Pi d% 
p 2 dK' 1 



K = K ffi + \Ze\ 



p 



2m„K + K 2 



18 Formally, it is given by 

oo 

R{K) = ]T Jo (C^) exp 



Konv-^ 

2K.(K) 



Vn{L) 



A n smh(S n L/2) 



with 



Vn{Z) 



Jo 



dr r Jo((^ n r / R) / dzexp 



V conv (Z — z) 
L 2/C(X) 

K,{K) S n coth(S , „i/2) with S, 



smh(S n (Z - z)/2) 



p(r, z) 



(31) 



(28) 



(29) 



,1/2 



encode 



The coefficients A n = 2hT ann +V com +K.(K) S n coth(S , »i/2) with S n = (V? onv /K,(Kf + Kl/R 2 , 
the effects of diffusion. 

19 The fit functions reproduce the results of our numerical calculation to better than 6% (with the excep- 
tion of the EinastoB MED case (for annihilations), for which the accuracy drops to a still acceptable 11%) 
over the whole range K = 100 MeV —> 100 TeV. The fit functions should not be used outside of this range. 
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Figure 9: Propagation function for antiprotons 
from annihilating DM, for the different halo profiles 
and sets of propagation parameters, and the corresponding 
fit parameters to be used in eq. (30) . 
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Figure 10: Propagation function for antiprotons 
from decaying DM, for the different halo profiles and 
sets of propagation parameters, and the corresponding fit 
parameters to be used in eq. (30). 
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The so-called Fisk potential (ftp parameterizes in this effective formalism the kinetic energy 
loss. A value of (ftp = 0.5 GV is characteristic of a minimum of the solar cyclic activity, 
corresponding to the period in which most of the observations have been done in the second 
half of the 90's and at the end of the years 2000's. 

4.1.5 Antideuterons 

The propagation of antideuterons through the Galaxy follows closely that of antiprotons 
discussed above, with a few trivial changes. The diffusion equation is still the one in eq. (24). 
In it: 

- Diffusion, being governed by the electromagnetic properties of the particles, is the 
same for antideuterons as for antiprotons, but of course the deuteron mass rrid should 
replace the proton mass in the expression for the kinetic energy and the momentum. 

- It is actually customary for low Z nuclei to use as a variable Kd/n: the kinetic energy 
per nucleon (n = 2 in this case). We will present all results as functions of this 
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Figure 11: Propagation function for antideuterons 
from annihilating DM, for the different halo profiles and 
sets of propagation parameters, and the corresponding fit 
parameters to be used in eq. (32). 
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Figure 12: Propagation function for antideuterons 
from decaying DM, for the different halo profiles and 
sets of propagation parameters, and the corresponding fit 
parameters to be used in eq. (32). 
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quantity. 

- The treatment of spallations of d on the interstellar gas ('annihilating' and 'non- 
annihilating' reactions) is less straightforward than for p, essentially for the scarcity of 
experimental nuclear data on d. We still write r( non _) ann = (nn + 4 2 / 3 nHeVpj° n axm v d 
and so we now need cr^ el . This can be obtained from related experimental measure- 
ments with the charge conjugated reaction pd or with the reaction pp: we refer for 
more details to [124, 125] and references therein. All in all, we find that a good 
approximation is to effectively adopt cr^j 1 ~ 2 cr^f 1 . 

With the ingredients above one can compute, exactly as for antiprotons, an antideuteron 
propagation function R d (K d /n) for annihilations and for decays for any choice of DM galac- 
tic profile and for any choice of set of propagation parameters among those in Table 1. We 
provide Rd{K d /n) for all these cases in terms of a fit function 

log 10 [R d (K d /n) /Myr] = a + at n d + a 2 K d + a 3 k\ + a 4 K d + a 5 K d , (32) 
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with Kd = log 10 (Kd/n/ GeV) and the coefficients reported in the tables in fig. 11 and 12 (and 
also reported on the website [29]). 20 Not suprisingly, since the changes are so minimal with 
respect to antiprotons and affecting subdominant processes only, the propagation functions 
resemble those for antiprotons closely. 

With these ingredients, it is straightforward to compute the antideuteron differential 
flux at the position of the Earth as 

1 dNl 
R d (K d /n) 2^ f dK~ ( anmmlatlon ) 

' dN( ■ < 33) 



dE 



47T 



PQ 

P& 
Mi 



DM 



\ dm 

\R d (K d /n)Y,T f ^ (decay) 



(notice that the primary d fluxes at production are given in Sec. 3 as a function of K d and 
not Kd/n). Solar modulation can be applied as discussed for antiprotons, with the due 
changes. 



4.2 Fluxes after propagation at the location of the Earth 

In this section we present the resulting fluxes of charged cosmic rays from DM as they 
would be observed at Earth, computed applying the formalism of the previous section. 



4.2.1 Electrons or positrons 

Applying the recipe of eq. (13) it is straightforward to compute the fluxes of electrons and 
positrons at Earth. We provide them in numerical form on the website [29], both in the 
form of Mathematica® interpolating functions and numerical tables. 

Fig. 13 presents some examples of such fluxes, for the cases of annihilation and decay. 
The choice of propagation parameters (MIN, MED or MAX) affects the final results up to 
one order of magnitude, especially at energies below tens of GeV and determine somewhat 
the spectral shape. For annihilating models, the choice of DM halo profile has a limited 
impact and again it gives rise to a spread in the prediction mainly at small E. Instead 
it is imperceptible for the decaying case. The black solid lines in the left panels of fig. 13 
represent the positron flux at Earth computed with the approximated method discussed in 
sec. 4.1.3 (and considering MED propagation parameters). As one can see, at the location 
of the Earth there are some differences in the flux up to a factor 2 due to the different shape 
of the halo function, but not much more than that. Therefore we can infer that the new 
full computation of the propagated flux has a limited impact, making the old one quite 
safe. Important differences, instead, are appreciable in other regions of the sky, especially 
close to the galactic center where one has larger energy losses. As a consequence we expect 
that the diffuse 7 rays, produced by these propagated electrons/positrons everywhere (see 
sec. 6), will be more sensitive to the difference between the two methods. 



4.2.2 Antiprotons 

Applying the recipe of eq. (27) it is straightforward to compute the fluxes of antiprotons 
at Earth, for a given choice of halo profile and propagation parameters. We provide them 

20 Thc fit functions reproduce the results of our numerical calculation to better than 5% (with the ex- 
ception of the EinastoB case, for which the accuracy drops to a still acceptable 10%) over the whole range 
Kd/n — 50 MeV — > 50 TeV. The fit functions should not be used outside of this range. 
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Figure 13: Fluxes of electrons or positrons at the Earth, after propagation, for the case 
of annihilations (top row) and decay (bottom row). In the left panels the propagation parameters 
are variated, while the halo profile is kept fixed. The opposite is done for the right panels. The 
choices of annihilation or decay channels and parameters are indicated. 



in numerical form on the website [29], both in the form of Mathematica® interpolating 
functions and numerical tables. 

Fig. 14 presents some examples of such fluxes, for the cases of annihilation and decay. 
We do not correct for any solar modulation. It is apparent that the choice of propagation 
parameters (MIN, MED or MAX) affects in a relevant way the final result, up to a couple 
of orders of magnitude, even if the spectral shapes are not sensibly modified. The choice 
of the DM halo profile, instead, has a limited impact and it is barely visible for the decay 
case. This is already evident of course in the little variations of the halo function in Fig. 10 
and can be traced back to the fact that the decay signal, being proportional to the first 
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Figure 14: Fluxes of antiprotons at the Earth, after propagation, for the case of an- 
nihilations (top row) and decay (bottom row). In the left panels the propagation parameters are 
variated, while the halo profile is kept fixed. The opposite is done for the right panels. The choices 
of annihilation or decay channels and parameters are indicated. 



power of the DM density, is mainly sensitive to the local DM halo, where the profiles do 
not differ sensibly 

4.2.3 Antideuterons 

Applying the recipe of eq. (33) it is straightforward to compute the fluxes of antideuterons 
at Earth, for a given choice of halo profile and propagation parameters. As usual, we 
provide them in numerical form on the website [29], both in the form of MATHEMATICA® 
interpolating functions and numerical tables, but we do not present example plots in this 
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case. The qualitative features are very similar to the case of antiprotons, with the necessary 
changes in scales. 



5 Prompt gamma rays 

Dark Matter produces high energy gamma rays both by direct ('prompt') emission during 
annihilation or decay and by Inverse Compton Scattering of produced by DM on the 
ambient light ('secondary'). In this section we focus on prompt gamma rays, while in the 
next one we deal with secondary emissions. 

The differential flux of photons from a given angular direction dfl produced by the 
annihilation of self-conjugated DM particles (e.g. Majorana fermions) is 

d$ 7 lr / p V T v^ ; . dNf f ds fp(r(s,e))\ 2 , 

— 3- = — ° j\ / av ).—JL J = — V V ;; annihilation 

dtldE 24n\M DM J ^ U dE 1 ./,._/•. V Pe J 

(34) 

where dN^/dE is the energy spectrum of photons produced per one annihilation 21 in the 
channel with final state /. If DM is not constituted by self-conjugated particles (e.g. in the 
case of Dirac fermions), then o~v must be averaged over DM particles and antiparticles: in 
practice, the equation above has to be divided by an additional factor of 2 if only particle- 
antiparticle annihilations are present. 

In the case of DM decay, an analogous equation holds 

MdE = T,M^ J ^ T ^ ' L.r.{ /, J ^ (35) 

Here the coordinate r, centered on the Galactic Center, reads r (s, 8) = (r^+s 2 — 2 r s cos d) 1 ^ 2 , 
and 9 is the aperture angle between the direction of the line of sight and the axis connecting 
the Earth to the Galactic Center. 



5.1 J factors 

The J factor in eq. (34) and eq. (35) integrates the intervening matter along the line of 
sight (along which the variable s runs) individuated by the angular direction, and it is 
conventionally weighted by r (here assumed to be 8.33 kpc) and the appropriate power of 
Pq (here assumed to be 0.3 GeV/cm 3 ) so to be adimensional. 22 J(8) is of course invariant 
under rotations around the axis which connects the Earth to the GC, due to the assumed 
spherical symmetry of the DM distribution p(r). 

The J factors are plotted in fig. 15 as a function of 9. We provide them in terms of 
M athematica® interpolating functions on the website [29] . 

The recipes (34) and (35) are ready for consumption if one needs the flux of gamma rays 
from a given direction. More often, of course, one needs the integrated flux over a region 
AQ, corresponding e.g. to the window of observation or the resolution of the telescope. 
The J factor is then replaced by the average J factor for such region, simply defined as 

21 Not per initial state particle; not per final state primary particle. 

22 Alternatively, sometimes an analogous factor is defined as J = p 2 (r) — r^pf^J in units of 

GeV 2 /cm 5 (annihilation) or J = f { p(r) = r 0/ 9 J in units of GeV/cm 2 (decay). 
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Figure 15: J(0) for annihilating (left) and decaying (right) Dark Matter, for the different 
DM profiles. The color code individuates the profiles (Burkert, Isothermal, Einasto, EinastoB, 
NFW, Moore from bottom to top in the inset). 



J (AO,) = (Jaq J dn) /An. The following simple formulae hold for regions that are disks of 
aperture # max centered around the GC, annuli 9 min < 9 < 9 max centered around the GC or 
generic regions defined in terms of galactic latitude b and longitude £ 23 (provided they are 
symmetric around the GC): 

dO sin 9 J (9), (disk) 

d9 sin 9 J (9), (annulus) 

f dbdi cosbJ(9(b,£)), (bx£ region) 

(36) 

where the integration limits in the formulae for J are left implicit for simplicity but obviously 
correspond to those in An. For the l b x I region' the limits of the integration region are 
intended to be in one quadrant (e.g. the b > 0°, < £ < 90° one for definiteness) , hence 
the factor of 4 to report it to the four quadrants. 

The values of the J factors and An for some popular observational regions are reported 
in table 2, for the cases of annihilating and decaying DM and for the different halo profiles. 
Any other region can be computed by using the formulae in eq. (36) and the J(9) functions 
provided above. 

23 Galactic polar coordinates (d,£,b) are defined as 

x = d cos i cos b, y = d sin I cos b, z = d sin b 

where the Earth is located at x = (such that d is the distance from us); the Galactic Center at x = r®, 
y = z = 0; and the Galactic plane corresponds to z w 0. Consequently cos0 = x/d = cos b ■ cos£. 
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With these ingredients, one explicitly has for the differential 7 ray flux from a region 



AQ 



dE 



Air 



( Pq V - ^ dN l 



Pe 



Mi 



DM 



f 

dNf 
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(decay) 



6 Photons from electrons and positrons in the Galaxy 

Galactic e generated by DM in the diffusion volume lose essentially all their energy into 
photons by means of two processes: Inverse Compton and synchrotron radiation. 

The resulting fluxes of ICS 7 rays and of microwave synchrotron radiation are thus 
possible signatures of DM. The ICS flux is particularly promising. One of its best features 
is that it originates from 'everywhere' in the diffusion volume of the galactic halo, including 
regions where the astrophysical background is reduced (e.g. at high latitudes). Moreover, 
essentially everywhere synchrotron energy losses are sub-dominant with respect to Inverse 
Compton energy losses (as discussed in Sec. 4. 1. 1), so that, thanks to energy conservation, 
the resulting ICS 7 flux suffers only moderate astrophysical uncertainties. 

The microwave synchrotron emission, on the other hand, is generated in significant 
amount from where the intensity of the magnetic field and of Dark Matter is highest, close 
to the Galactic Center, and therefore is plagued by more uncertainty and more background. 

We first present the detailed recipe to compute ICS 7 rays. We leave synchrotron 
radiation for the next subsection. 



6.1 Inverse Compton gamma rays 

The differential flux of ICS photons within an angular region AQ can be written in terms 
of the emissivity j(E^, r) of a cell located at a distance r = \x\ from the Galactic Center as 

dE 1 dVL £ 7 y LoA 4vr v ; 

In general, for any radiative process, the emissivity is obtained as a convolution of the 
spatial density of the emitting medium with the power that it radiates (see e.g. [126]). In 
this case therefore 

rM DM (/2) dn 
j(E^r) = 2 dE e P IC (£ 7 ,£ e ,r) —^( r ,E e ), (39) 

where V = V\ c is the differential power emitted into photons due to ICS radiative 
processes (the sum runs over the different components of the photon bath: CMB, dust- 
rescattered light and starlight) and dn e ±/dE e is the electron (or positron) number density 
after diffusion and energy losses, as computed in subsection 4.1.2 (notice that there it was 
denoted as / for simplicity, see page 18; dn e ±/dE e just corresponds to eq. (13) removing 
the v e ± 1 47r factor) . The minimal and maximal energies of the electrons are determined by 
the electron mass m e and the mass of the DM particle M DM - The '/2' notation applies to 
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the decay case. The overall factor of 2 takes into account the fact that, beside the electrons, 
an equal population of positrons is produced by DM annihilations/decays and radiates. 24 

The radiated power Pic, in the full Klein-Kishina case, is given by (we refer the reader 
to [106, 107] and references therein for more details on the derivation) 



3^tA/, E, \n i {E°(q),x) 

'l/4 7 2 



47 2 A/^r 7 4g 7 2(l- e) ; q 



o 1 e 2 1 (40) 

2glng + q + 1 - 2q 2 + ~r^(l - q) 



where 7 = E e /m e is the Lorentz factor of the scattering electron and the integrand is 
expressed in terms of 



4£°£ e _EL L . J. 
r E (l-e)' ' ml ' C " E e ' m 4 7 2 



with T E = ^ 7 9 e , e = in ~ < g < 1. (41) 



Here E® is the initial energy of the photon in the background bath. Correspondingly, E 1 
lies in the range E®/E e < E 1 < E e YE/{\ + T^). The non-relativistic (Thompson) limit 
corresponds to < 1, so that e < 1, the last term in the integrand of P is negligible, 
and q — > y = E 1 /{A^ 2 E^ j ) with < y < 1. Thus in the Thomson limit 

V l lc (E^ E e ,x) = ^E 1 [ dy Ui ( E ^> x } [ 2 y In y + y + 1 - 2y 2 ] [Thomson limit]. 

47 Jo y 

(42) 

Plugging now Pic and n e ± in eq. (39), we can write the IC differential flux in the 
following convenient form: 



d^ic-y 1 r Q 
dE 1 dVl ~ £ 2 Att 



2 VM DM 
M DM 



1 1 / o \ 2 f M ° M dN^ 

M dE s J2(™)f^(E s )I lc (E^E s ,bJ) (annihilation) 

M DM /2 ^jy/ 

r / ^(E s ) /ic^, £ s , 6, t) (decay) 

(43) 

where E B is the e 1 * 1 injection energy and I\c{E 1 , E s , b, £) (with the dimension of an energy) is 
a halo function for the IC radiative process. This formalism allows therefore to express the 
flux of ICS 7 as the convolution of the electron injection spectrum dN e ±/dE and this new 
kind of halo functions, in close analogy with the formalism for charged particles. Indeed, we 
can explicitly express lie i n terms of the ICS ingredients discussed above and the generalized 
halo functions for e ± that we introduced in Sec. 4.1.2. We get 



lies. r Q \ P& J Jm e K E , r ( s , & )) 

(44) 

where again rj = 1, 2 for the decay or annihilation scenarios respectively. The intensity 
of the interstellar radiation Y2i n i cancels out in the ratio ^2Vi/b, up to the sub-leading 
synchrotron contribution and provided that we are not interested in the contributions from 
the individual light baths. 



24 Recall from footnote 9 that with the notation we always refer to the independent fluxes of electrons 
e~ or positrons e + and not to the sum. 
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The class of functions Iic(E 7 , E s , b, £) are plotted in fig. 16 for the annihilation case (and 
in fig. 17 for the decay case) by varying the e ± injection energy and the latitude b for a fixed 
longitude £ = 0.1°. As one can see, increasing the latitude (from 0.1° to 90°) the spectrum 
becomes less step due to two combined effects: the decreasing abundance of star-light and 
dust, that provide the more energetic part of the target photon bath, and the always present 
Klein-Nishina relativistic effect, that blunts the high energy part of the spectrum especially 
for large injection energy and in regions close to the GC where more star-light and dust 
are present. Furthermore, at small latitudes, the difference in the normalization among 
different profiles becomes appreciable. 

All these classes of functions for decay and annihilation scenarios are provided in 
M athematica® interpolated functions form at the website [29] . There is one halo function 
for each DM density profile and e propagation model (MIN, MED, MAX); each function 
depends on 4 variables: the energy E® of as produced by DM before diffusion; the final 
photon energy Ey; the two galactic coordinates b, £ that specify the observed direction. 

The final step to obtain the differential ICS 7 flux d$ic 7 / 'dE^dVt consists in performing 
the convolution integral over E s with any desired prompt energy spectrum from DM. 
We do not provide pre-compiled ICS 7 fluxes (due to file-size limitations), but we provide 
the code that performs such final integral on the website [29]. 

Finally one can compute the differential flux from a region AQ by integrating over b 
and £ as discussed in Sec. 5 



Due to the intertwined dependence on b and £ and on E 7 and E s of the halo functions, 
here the geometrical integral cannot be pulled out as for prompt 7 rays, so a J factor cannot 
be defined in a simple way. 

Figure 18 presents some examples of IC 7 ray fluxes computed as discussed above. In 
the top left panel (for the case of annihilations into bb) we plot the integrated flux in a 
region which is large but quite close to the GC, and we vary the propagation parameters 
of the electrons and positrons: the impact on the IC 7 fluxes is small but visible at small 
energies, and it is due to the difference in the resulting populations of e ± . In the top right 
panel we vary the DM profiles in the '5x5' region: since the profiles differ sensibly in such 
a relatively small region, the impact on the fluxes is of the order of more than an order of 
magnitude. The lower panels refer to decay cases and focus on leptonic channels. When 
looking for the signal in the region of the Galactic Poles, the variation of the propagation 
parameters (left plot) has a significative impact mainly because the height of the diffusion 
cylinder (and therefore the contribution to the integrated emission) increases as one moves 
from MIN to MAX. When looking for the signal from a small '3x3 region around the GC 
(right panel), the variation of the choice of profile has a non-negligible impact (even in the 
decay case) because the observation is concentrated towards where profiles differ most. 

6.2 Synchroton radiation from e ± 

In this section we briefly discuss the synchrotron emission from DM produced energetic 
e immersed in strong magnetic fields (see e.g. [127] and references therein). We do not 
provide, however, numerical results for this signature. 




(45) 
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Figure 16: Halo functions for Inverse Compton Scattering, for the Dark Matter annihi- 
lation case and varying DM halo profiles and electron injection energies. 
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Figure 17: Halo functions for Inverse Compton Scattering, for the Dark Matter decay 
case. 
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Figure 18: Fluxes of galactic Inverse Compton gamma rays, for the case of annihilations 
(upper panels) and decay (lower panels). 



An electron or positron with momentum p in a turbulent magnetic field B generates 
the following spectrum of synchrotron power W syn into photons with frequency v: 



dW syn V3 e 3 B 



dv 6ir m 



F(— ), F(x) = x K 6/i {Z)d£ « ^=5(x - 1/3) (46) 



so that, to a good approximation, it generates photons with frequency u syn /3, where 

3eBp 2 B ( p x2 



^=4^F= 4 - 2MHZ G^J <47) 

such that the total energy radiated is proportional to B 2 . Galactic magnetic fields are highly 
uncertain; we adopted the simplest 'conventional' choice of eq. (10). While we included in 
eq. (8) the magnetic field contribution to e energy losses, it is sub-dominant everywhere 
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with respect to the Inverse Compton contribution. Consequently, as long as one considers 
any other magnetic field profile which is similarly sub-dominant, it is not necessary to 
recompute in each case the halo functions describing the diffusion of e ± . (The diffusion 
coefficient, which also depends on the turbulence spectrum magnetic field, is in practice 
extracted by fitting models of cosmic rays as discussed above). 
Consequently a galactic distribution dn e ± / dE e generates 



dW f 



syn 



dv dVt 



2 

47T 



ds 



l.o.s. 



e dE P dv 



(4f 



Inserting dn e ±/dE e from (13) and dW syn /dv from eq.s (46) and using the ^-function ap- 



proximation such that E e 



V 



yjAirnilv/eB = 0.43 GeV(z//GHz) 1 / 2 ( J B/mG)- 1 / 2 we get 
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where e is the electron charge and J-^ = ub/u is the fraction of the e ± energy lost into 
synchrotron, with an uncertain spatial profile (u is defined in eq. (8)). 



7 Extragalactic gamma rays 

The 7-rays emitted by DM annihilations or decays in all the extragalactic structures and 
(in principle) all along the history of the Universe reach us in the form of an isotropic 
contribution to the total 7-ray intensity. In this section we discuss how to compute it, 
in terms of many of the ingredients introduced above. With respect to the galactic case, 
however, we have to include the effect of the 'cosmological dimming' due to the expansion 
of the Universe and the fact that, unlike in the galactic environment, on cosmologically 
large distances one can not neglect the absorption of gamma-rays. We mainly follow the 
formalism presented in [128, 68] (see also [129]). 

In full generality, the differential flux of extragalactic gamma rays perceived at a certain 
redshift z is given by 

00 3 

l"wh*) G£) wb; ' 2,) e " I( ""''- (50) 

z 

The factor [H(z')(l + z')]" 1 , where H(z) = H y/n m (l + zf + (1 - VL m ) is the Hubble 
function, converts the redshift interval to a proper distance interval (in other words, the 
integral can be thought to be the integration over time of all the photons emitted in the 
past, then converted into redshift). The factor [(1+z)/ (1+z')} 3 accounts for the cosmological 
dimming of intensities due to the dilution of the source. E' is the energy of a photon at 
redshift z', assuming it has energy E at redshift z: E' = E(l + z')/(l + z). The function 
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t(E^, z, z') is the optical depth describing the absorption between the redshifts z and z' . 
To calculate the Hubble function we assume a flat ACDM "concordance" cosmology with 
tt m = 0.27, tt A = 0.73 and h = 0.7. 

Because we observe the gamma rays at z = 0, the formula for the flux simply reduces 

to 
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dE*, 
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dz' 
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H(z')(l + z')' 



-t{E 1 ,z') 



(51) 



As in the Galaxy, the local gamma-ray emissivity for the annihilating/decaying DM 
models is the sum of the prompt and IC radiation parts 
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The emissivity jvG™ pt is given by 
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(53) 



where the average cosmological DM density is p(z) = po(l+z) 3 and po — 115 10~ 6 GeVcm -3 
is its value today [z = 0); dN 1 /dE 1 is the spectrum of the prompt photons as computed in 
section 3. In the case of annihilations, DM clustering enhances the rate: the factor B(z) 
effectively takes that into account (as described in detail below). 

The computation of the emission coefficient for Inverse Compton radiation is more 
involved: IC radiation depends on the population of the e ± produced by annihilation/decay 
of DM at redshift z' and the background photon bath n(E', z') at the same redshift. Ac- 
tually, more precisely, the population of the e ± at a certain redshift z' is the result of the 
diffusion-and-loss processes (in space and 'time') from preceding redshifts. However, the 
mean free path of in the intergalactic medium is very short compared to cosmological 
length scales (dominantly due to the presence of the CMB). Thus we can approximate the 
IC spectrum as generated "on-spot" and neglect the effect of such diffusion. (Technically, 
remember that this corresponds to setting the equivalent of the 'generalized halo function' 
defined in Sec. 4.1.2 to 1.) In turn, concerning the background photons, it is a good ap- 
proximation to assume that they only consist in the CMB: one can safely neglect the IR, 
light and UV photons from starlight and secondary dust radiation that were instead impor- 
tant in the Galactic medium, where ICS on CMB, IR and light photons give comparable 
contributions. 

Under these approximations the IC emissivity can be expressed as 
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where the overall factor of 2 takes into account that equal populations of electrons and 
positrons radiate (the '/2' notation applies to decay). The functions V^ AB (E' , E e , z') and 
bfo [B (E e , z') are the radiated power and the energy loss coefficient function for e , exactly 
analogous to those defined in Sec. 6.1 (eq. (40)) and Sec. 4.1.1 (eq. (9)), the only difference 
being that they are now functions of the redshift and not of the position in the Galaxy. They 
can be computed in the full Klein-Nishina case or in Thomson limit. As the intergalactic 
medium is dominated by low energy CMB photons, the Klein-Nishina formalism is needed 
only for the extreme mass region of DM, above Mdm > 20 TeV (see the comparison between 
the black and colored dotted lines in the last panel in Fig. 20) (the situation is therefore 
different from the case of the Galaxy, where IR and stellar light gives a relevant contribution 
to the target photon bath, and thus the Klein-Nishina formalism should be used already at 
energies above M DM ~ 100 GeV). 

The advantage of using the Thomson approximation, when valid, is that one needs 
to calculate the IC emissivity only for one fixed redshift, due to the fact that n(E, z) = 
^cmb(-^, z) = n Pi[E, T (l + z)), where n P1 is the black body spectrum and T is the tem- 
perature of CMB at z = 0. For example, one can calculate the IC spectrum at z = and 
then for other redshifts use simple scaling properties of n-p\. Heuristically, this can just be 
understood as follows: the IC cross section in the Thomson regime is energy-independent; 
since the energy spectrum of the injected e is the same at any z, the IC photon at the 
production epoch z will have an energy proportional to the one of the upscattered back- 
ground photon. The latter is (1 + z) times the current one, which is exactly the factor 
compensated by the subsequent redshifting. So, the IC spectrum is universal and equal 
to the one calculated at z — 0. In formulae, this means that the ICS emissivity is just 
given by j^~(E' z') = B(z')(l + z') 6 j^L (E 7 , 0) for the case of annihilating DM and by 
j^ 7 (E;, z') = (1 + z'f ]f Gl (£ 7 , 0) for the case of decaying DM. 

When the Thomson approximation is invalid and one has to apply the full Klein-Nishina 
formalism for the calculation of IC radiation, the emission coefficient of IC radiation J IC 
has to be calculated for all the relevant region of redshift. 

7.1 Effect of DM clustering 

The emission coefficients for annihilating DM described above contain a cosmological boost 
factor B(z) which multiplies the contribution from the homogeneous distribution of DM 
Pdm(z)- Technically, the true square density that enters the formulae is thus (pdmW) = 
p 2 DM (z)((l + S(z)) 2 ) = B(z)p 2 DM (z). Here 5 = p/p — 1 is the density contrast and the 
cosmological boost factor B(z) = ((1 + S(z)) 2 ) = 1 + (5 2 (z)). To calculate the boost factor 
we adopt a halo model which approximates the matter distribution in the Universe as a 
superposition of DM halos (e.g. [130]). Within this model B(z) can be given as 



where p m $ is the matter density at z = 0, A c ~ 200 is the overdensity at which the halos 
are defined and M m i n is the minimum halo mass. ^(M, z) is the halo mass function, that 
can be cast in the universal form [131] 



B(z, M min ) = 1 + 
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where the parameter v = [5 c (z)/a(M)] 2 is defined as the ratio between the critical overden- 
sity S c (z) and the quantity a(M) which is the variance of the linear density field in spheres 
containing a mean mass M. For the multiplicity function f{y) we use the form in [132] 

i 

vf{v) = A (l + -1 )(f\* e-"P (57) 



where v' = av with a = 0.707, p = 0.3 and A is determined by requiring that the integral 
J dvjiv) = 1. Eq. (57) reduces to the original Press-Schechter formula [131] taking a — 1, 
p = and A = 1/2. 

c(M, z) represents the halo concentration parameter function and the function /(c) for 
the halos with the NFW density profile [133] is given as 



c 3 



1 



1 + c 



log(l + c) 



1 + c 



(5* 



For the concentration parameter function c(M, z) we use two different models: 



• the 'Maccid et aV model [134], in which c[M,z) = k 200 (U{z f (M)) /H(z)) 2/3 , where 
&200 — 3.9, H(z) = H(z)/H and Zf(M) is the effective redshift for the formation of 
a halo with mass M; 

• a 'power law 7 model (inferred from the results in [135, 134]) with c(M, z) = 6.5 H(z)~ 2 / 3 
(M/M*) -ai , M* = 3.37 1O 12 /i _1 M , which gives a good fit within the mass range re- 
solved by the simulations. 

Moreover, we consider two typical values for the minumum halo mass (motivated e.g. 
by [136, 137]) : 

• M min = 1O~ 6 M and 

• M min = 1O- 9 M . 

The resulting four models for the cosmological boost factor are given in form of a Mathema- 
tica® interpolating function on the website [29] and are plotted in Fig. 19a. We see that 
the boost factor shoots up at z ~ 100, which corresponds to the redshift where the DM 
halos with the assumed lowest masses start to form. As expected, the curves for the 10 -9 
M & models start to deviate slightly earlier. 

It is clear from the figure that currently the biggest uncertainty influencing the mag- 
nitude of the boost factor is related to the concentration parameter function c(M, z) and 
in particular to the way one chooses to extend its behavior below the mass scales directly 
probed by the simulations. At low redshifts the variation between the models reaches two 
orders of magnitude. On top of that there is an extra uncertainty related to the profiles 
of the DM halos: we have assumed the NFW profile as a benchmark, but expressions 
analogous to eq. (58) (but more cumbersome) can be easily calculated for the Einasto and 
Burkert profiles. Comparing the f(c) functions for different profiles along a large range of 
Cvir, we find that this adds an uncertainty of roughly one order of magnitude. 
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Figure 19: Left: Extragalactic (cosmological) boost factors B(z) for various concentration 
models and for various lower cut-offs of DM halo mass. Right: Contour plot of the opacity of the 
Universe, for the 'min UV background case (color plot). The short-dashed curves correspond to 
'no UV background case. 



7.2 Absorption of gamma rays in the intergalactic medium 

The factor exp[— r(E 1 , z, z')] in eq. (50) expresses the attenuation of the flux of 7-rays 
having an energy at collection E 7 = E' (1 + z)/(l + z') as they propagate from the emission 
redshift z' (where they are emitted with energy E') to the collection redshift z. In turn, the 
notation r(E 1 , z') identifies the optical depth for gamma rays collected today with energy 
Ej = E' /{I + z'). In this section we briefly discuss the relevant physics processes and 
provide the ingredients to compute t(E 7 , z, z') explicitly. We provide t(E 7 , z') in the form 
of a Mathematica® interpolating function on the website [29]. 

The processes relevant to the absorption of energetic photons in cosmological length 
scales (> 10 Mpc) and in the energy range roughly spanning from MeVs to TeVs are 

- pair production on baryonic matter 

- photon-photon scattering on ambient Photon Background Radiation (PBR) 

- pair production on ambient PBR 

Fig. 19b illustrates the value of e~ T as a function of the energy at detection E^ and the 
emission redshift z' . In the lowest energy section of the plot, before the first small plateau 
located at KeV-ish energies, the absorption is dominated by photoionizations. Beyond that, 
up to the beginning of the large plateau, Compton losses are dominating. In the large flat 
plateau region the absorption is dominated by pair production on matter. In the final falling 
part of the curves the absorption is determined by photon-photon pair production. The 
photon-photon scattering also gives some contribution in the region where the flat plateau 
turns over to falling curves. 
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The PBR is mainly composed by the CMB, the intergalactic stellar light and secondary 
IR radiation. The intergalactic stellar light notably consists of the UV background produced 
in the low redshift Universe once the first (massive and hot) stars start to light up. Since 
this latter part is the most uncertain, we introduce three distinct modelizations: 

• 'no UV assumes that no UV background is present. 

• 'minimal UV takes into account that recent studies of blazars, e.g. [138], suggest 
signicantly lower values for the UV photon densities than estimated in many of the 
previous investigations. We use the UV background model as given in [139], which is 
fully consistent with the earlier study [140]. 

• 'maximal UV assumes the UV background as given by 'minimal UV is multiplied 
with factor 1.5. It can considered as an upper limit of possible systematic errors of 
the model given by [139]. 

We will see that the choice of UV background has a certain impact on the flux of very 
high energy extragalactic gamma rays. 

We now proceed to discuss a compact formalism allowing to calculate the optical depth 
t{E' z,z), focussing in particular on the last three processes in the list at page 46, since, as 
discussed, they are most relevant for GeV gamma ray astronomy telescopes such as FERMI. 
For further details see [141] and [142]. In what follows we express the photon energy in 
units of electron rest mass, i.e. e = i? 7 /(m e c 2 ). 

o The absorption coefficient at redshift z for pair production on neutral matter with 
mass fractions of hydrogen X = 0.75 and helium Y = 0.25 can be approximated as 



Here f2& = 0.0448 ± 0.0011 is the density parameter for baryons [7], h is the reduced 
Hubble parameter, «/ = 7.29735 x 10~ 3 is the fine structure constant, and r = 
2.8179 x 10~ 13 cm the classical electron radius. 

Below redshift z ~ 6 the Universe is reionized [143]. The absorption coefficient for 
pair production on fully ionized matter can be approximately given by 




(59) 



(60) 



where the average electron number density at z 



is 




(61) 




a (l + z) 3 ln(2e) 
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(62) 
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e > 1 , 




(63) 
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o The absorption coefficient for photon-photon scattering at redshift z can be given by 

^77— scat ( 



a 7 7-scat(e, z) = a (l + z) 6 e 3 , (64) 



= 4448^ / I^y 

455625 A c V 2 -725Ky P ' K ' 



where Go = is the CMB temperature at z = in electron rest mass units, 

and A c = 2.4263 x 10 -10 cm is the electron Compton wavelength. Here we consider 
the CMB only as target photons: the contribution of the UV photons from stars is 
completely negligible. 

o For photon-photon pair production processes, instead, the contribution of the inter- 
galactic UV background is not negligible, so we have to include those together with 
the CMB. The target photon number density is thus n(e, z) = ucmb^, z) + n V y(e, z). 
As there are no simple analytical forms available for nuv(e, z), we have to use a fully 
numerical treatment for the photon-photon pair production process. The absorption 
coefficient for photon-photon pair production at redshift z is then computed as 

00 



X 2 c a 2 f f <f)(v) 
a 7 7- P air(e, z) = — — / den(e, z)— ^, (66) 



l + 2f + 2f 2 1 2^(1 + 2^) , 2 , 2 , 
with 6(v) = Inw-^^ '- - In 2 w + 21n 2 (l + w) + 
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+ 4 ^{rh)-T- (67) 

where v = (e e — 1) > , w = — V C . (68) 

V 1 + v — Jv 



Here Li2 is the dilogarithm function. 
With the above ingredients, the total absorption coefficient a(e, z) is given by 

a ma t-pair(e, z) if 6 < z < 1000 \ , s , , (m 

«ion ma t- P ai r (e, z) if z < 6 J + "77-scatle, *j + a 77 - pai rle, *J. (bJj 

So finally the optical depth r(_E 7 , 2, 2/) between redshifts z and z' can be computed as 



z 



„ l + z E~, 
where e 



l + z m e c 2 



7.3 Extragalactic gamma rays: results 

Applying the recipe of eq. (50), with all the ingredients discussed in this section, we com- 
pute the fluxes of extragalactic gamma rays detected at Earth, for any given choice of the 
concentration parameter function, minimal halo mass and UV background. We provide 
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Figure 20: Fluxes of extragalactic gamma rays, for the case of annihilations (first 3 panels) 
and decay (last panel). In each panel one of the astrophysical model assumptions is variated. The 
choices of annihilation or decay channels and particle physics parameters are indicated. On the 
last panel, the black dotted line indicates the flux for the 'no UV case computed in Thomson 
approximation. 



them in numerical form on the website [29] in the form of MATHEMATICA® interpolating 
functions. 

Fig. 20 presents some examples of such fluxes, for the case of annihilation and decay and 
for different primary channels. In the top left and top right panels the minimal halo mass 
M min and the choice of c(M, z) are variated respectively. As anticipated in the discussion 
of the cosmological boost factor (see 7.1), this affects the normalization of the spectra by a 
factor of a few (varying M m i n ) or almost by two orders of magnitude (varying the c(M, z) 
model). In the bottom panels the UV background (relevant in the gamma rays absorption, 
see 7.2) is varied. The effect is visible at energies above a few TeV. In the last panel we 
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also plot for illustraton the flux (for the 'noUV case) that one would obtain if adopting 
the Thomson approximation discussed in Sec. 7. For smaller DM masses the effect would 
be even less important. 



8 Summary 
8.1 Ingredients 

On the website [29], we provide the following: 

|TJ dlNdlxIEW[primary->f inal] [DMmass , log 10 x] : spectrum d In N/d In x in x = E/M^u 
of final particles generated by DM annihilations into a pair of primary particles. 
Also in terms of numerical tables. 



2j dlNdlxI [primary->f inal] [DMmass , logiox] : same as but without EW correc- 
tions. 

Also in terms of numerical tables. 



3 b[E,r,z]: energy loss coefficient function b(E,x) for e ± of energy E at the position 
(r,z) in the Galaxy. 



4j ElectronHaloFunctGalaxyAnnI [halo , propag] [log^x, logioE s ,r ,z] : generalized halo 
functions I(E, E s , r, z) for e for annihilations, in any given point (r, z) of the Galaxy, 
expressed function of x = E/E s . 
Analogous for decay. 



5j ElectronHaloFunctEarthAnnI [halo ,propag] [log^x, logioE s ] : generalized halo func- 
tions I{E,E sl r & ) for e for annihilations, at the location of the Earth, expressed as 
a function of x = E/E s . 
Analogous for decay. 



6 Tables of fit coefficients for the reduced halo functions for e for annihilations 
X(A,r ) at the location of the Earth. 
Analogous for decay. 

|~7] Tables of fit coefficients for the propagation functions for p for annihilations 
R{K) at the location of the Earth. 
Analogous for decay. 



8j Tables of fit coefficients for the propagation functions for d for annihilations 
R(Kd/n) at the location of the Earth. 
Analogous for decay. 



9j ElectronFluxAnn [primary, halo, propag] [DMmass, o"v,log 10 E e ] : differential flux 

dE 

at Earth. 

Analogous for decay. 

d§t 



10 | ProtonFluxAnn [primary , halo , propag] [mass,o"v,log 10 K] : differential flux — — ^ at 

dK 

Earth. 

Analogous for decay. 
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1 1 DeuteronFluxAnn [primary , halo , propag] [mass , av , log 10 Krf] : differential ffux — — 

dK d 

at Earth. 

Analogous for decay. 



12 | Jave [halo] [Logio#] : factor J{6) for prompt gamma rays. 
Analogous for decay. 



13 | I ICAnnI [halo , propag] [Log 10 i? s ,Log 10 E 7 , Log 10 £,Log 10 6] : halo functions 7ic(-E 7 , E s , b 
for Inverse Compton, for annihilations. 
Analogous for decay. 
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Code bite to compute 



ic 7 



dE*y 



15 I BoostF [z ,minhalomass , cmodel] : cosmological boost factor B{z) as a function of 

redshift z for a choice of M m i n and c(M). 

16 ETautEjZ' ,UVmodel] : optical depth of the Universe e T ( E ' z ^ for a choice of UV back- 
ground. 

17 EGgammaFluxAnn [primary , minhalomass , cmodel , UVmodel] [mass , av , IE] : differen- 
tial flux — ^° 7 . 

dE 1 

Analogous for decay. 

The MATHEMATICA® InterpolationFunctions and the code bite provided in 14 



have been produced with MATHEMATICA® version 7.0.0. The InterpolationFunctions 
are expected to be compatible with version 2 and any later version. The code bite employs 
solution methods included in Mathematica® version 4 and later. 

Table 3 lists the discrete choices for the variables employed in the functions above, 
together with a reference to the corresponding discussion in the text (when available). 

8.2 Recipes 

The main recipes for computing DM indirect detection signals are: 

I. Computing 6 : the differential flux of e ± at Earth: 

dE 



eq. (13) with _3_|, [l] and |_5_|. Or use [9_ 



II. Computing 6 with approximated energy losses: 

dE 

eq. (22) with [T] and [6 . 

III. Computing the differential flux of antiprotons at Earth: 

dK 



eq. (27) with [T] and \j] in eq. (30). Or use [W. 
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Variable 


Values or range 


Refer to 




eL, eR, /iL, /iR, rL, rR, 




primary 


q, c, b, t, 7, g,WL, WT, ZL, ZT, 

ill -1 C , llioc, ill 7(i , llonn . tiffin « llAnn « ±lc;nr^ « 

1 / r\ 1/11 1 /T \/ !v a \7 ^ // \T ^ T~ 

//fci ; ^/-t-} £/ / j v r tJ ; V r v r / 


Eq. (2) 


~F "1 Tl a ~\ 


P T"\ 0/ H 7 /P 1/11 1 /1~ 


Oct. o 


DMmass 


^ PpV —V 1 fin TpV fnrmihilntinn^ 

10 QpV — i 900 TpV Mpppiv^ 


Sec. 3.3 


11CLX VJ 


WFU Fin FiR Tcin Rnr Mnn 


Ficr 1 




MIN, MED, MAX 


Table 1 


crv or r 


any 




minhalomass 


1CT 6 , 1(T 9 


Sec. 7.1 


cmodel 


maccio, powerlaw 


Sec. 7.1 


UVmodel 


noUV, minUV, maxUV 


Sec. 7.2 



Table 3: Variables of the numerical functions and their admitted values. 



IV. Computing 



dK d 



the differential flux of antideuterons at Earth 



eq. (33) with \T\ and _8_ in eq. (32). Or use 11 . 

V. Computing -y=p-: the differential flux of prompt 7 rays: 

dE^ 



eq. (37) with \Y\ and J from table 2 or from eq. (36) with 12 



VI. Computing 



'IC7 



dE^ 



: the differential flux of galactic ICS 7 rays: 



eq. (45) with eq. (43), with \Y\ and 13 . Or use directly 14 



VII. Computing v- 



dv dVt 

eq. (49), with T and [4 



: the differential flux of galactic synchrotron radiation: 



VIII. Computing ^, G7 : the differential flux of extragalactic 7 rays: 

eq. (51) with eq. (53) and eq. (54), and with 15 , 16 and JL_. Or use directly 17 
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